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INTRODUCTION 


Essential  features  of  any  mathematical  method  can  be  studied 
best  if  extreme  cases  are  considered.  As  far  as  interpolation,  fil¬ 
tering,  prediction,  or  even  least-squares  collocation  are  concerned, 
an  infinite  homogeneous  set  of  regularly  distributed  data  presents 
itself  as  an  ideal  candidate  for  such  kind  of  studies.  Without  re¬ 
striction  of  generality  one  can  assume  a  unit  distance  between 
neighboring  data  points,  such  that  the  data  are  located  at  the 
places  of  the  cardinal  numbers  on  the  real  line;  this  is  why  we  have 
used  the  term  "cardinal"  in  the  title.  The  second  term  "interpola¬ 
tion"  does  not  require  a  further  interpretation,  although  the  pre¬ 
sent  study  goes  beyond  that. 

The  goal  is  to  investigate  the  response  of  the  choice  of  var¬ 
ious  base  (covariance)  functions  onto  the  interpolated  (predicted) 
function  based  on  an  infinite  homogeneous  data  set. 

Chapter  1  deals  with  the  whole  family  of  spline  base  func¬ 
tions,  starting  with  the  zero  degree  spline  base  function  (=  step 
function)  and  finally  arriving  at  the  highest  degree  spline,  whose 
corresponding  sampling  function  equals  the  sirntx/nx  function. 

In  the  second  chapter,  relations  between  Gaussian  functions 
of  various  correlation  lengths  and  the  corresponding  spline  func¬ 
tions  of  certain  degree  are  established.  Deviations  of  these  two 
functions  from  each  other  are  estimated. 

Structural  similarities  between  spline  interpolation  and 
least-squares  interpolation  are  discussed  in  chapter  3  .  The  fini¬ 
te  support  of  splines  and  their  ability  to  fit  and  replace  any- 
kind  of  covariance  (=  base)  function  leads  to  the  conclusion  that 
splines  of  high  and  odd  degree  may  advantageously  be  used  as  base 
functions  for  the  solution  of  least-squares  prediction  problems. 
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A  virtually  completely  different  interpolation  method, 
which  has  been  advocated  by  B^erhammar  and  others,  is  investi¬ 
gated  in  chapter  4  .  An  interesting  and  remarkable  relation 
to  spline  functions  of  odd  degree  is  presented. 

The  impact  of  data  noise  and  correlation  length  onto  the 
sampling  functions  and  the  prediction  error  is  investigated  in 
chapter  5  .  Of  particular  importance,  for  all  considerations 
concerning  data  collection,  is  a  very  strong  relation  between 
the  correlation  length  of  the  underlying  covariance  function, 
the  sampling  rate  (distance  between  data  points) ,  the  data  noise, 
and  the  maximum  prediction  error.  Its  impact  on  practical  gravity 
field  sampling  problems  might  be  considerable. 

There  exists  a  whole  bunch  of  base  functions  which  can 
be  used  for  interpolation  and  prediction  problems.  It  is  asto¬ 
nishing  and  remarkable  that  such  strong  ties  exist  between  all 
of  them.  The  beautiful  interplay  between  the  physical  and  the 
frequency  domain  is  given  a  dominant  role  in  all  investigations 
performed  here. 
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1 .  SPLINE  INTERPOLATION  ON  THE  REAL  LINE 


Let  us  give  a  definition  of  the  key  words  appearing  in 
the  above  title.  There  exist  basically  two  different  types  of 
interpolation  functions:  the  ones  are  based  on  the  intrinsic 
statistical  information  contained  in  a  data  set,  the  others  are 
not;  the  latter  ones  are  referred  to  as  "deterministic"  inter¬ 
polation  functions.  By  the  term  "local"  we  understand  an  inter- 
polant  which  is  sensitive  with  respect  to  data  in  the  neigh¬ 
borhood  of  the  interpolation  point  and  blind  or  considerably 
less  sensitive  with  respect  to  remote  data.  (A  typical  counter¬ 
example  would  be  a  single  interpolating  polynomial.)  An  interpo- 
lant  is  understood  as  a  function  which  reproduces  the  data.  When 
talking  about  "data"  we  have  a  homogeneous  set  of  error-free  data 
in  mind  which  is  regularly  distributed  along  a  line;  without  re¬ 
striction  of  generality  we  can  assume  the  data  to  be  located  at 
the  cardinal  numbers  .  ..,  -2,  -1,  0,  1,  2,  ...  ;  the  data  set 
is  supposed  to  have  infinite  dimension,  it  is  denoted  by  the  in¬ 
finite  vector  pair 

x  =  {j}  ,  f  =  { f  }  ,  j  =  ...,  -2,  -1,  0,  1,  2,  ... 

J  t 

with  x  denoting  the  vector  of  coordinates  and  f  the  vector 
of  corresponding  function  values. 

In  the  sequel  interpolants  will  be  studied  which  differ  by 
their  interpolation  and  differentiation  properties,  or  simply,  by 
rheir  smoothness. 


THE  STEP  FUNCTION  "INTERPOLATION 


In  many  numerical  integration  problems  che  most  primitive 
functional  representation  of  the  integrand  is  chosen,  essentially 
because  of  simplicity  (e.g.  Stokes'  integral);  it  is  the  step- 
function  representation  of  the  "true"  function  which  is  here  sam¬ 
pled  at  the  cardinal  numbers .  The  base  function  is  the  discontin¬ 
uous  and  strictly  local  function  (Fig.  1.1) 


U 

1 

for 

i  x !  <  2 

ll 

for 

X 

ii 

(1.1) 

!° 

else 

• 

CfiflC I f!fll  3  -  SPLINE  CF  DEGREE  0 


FIG  .  1.1  Base  function  for  K~ 


It  is  an  element  out  of  KQ  ,  the  space  of  discontinuous  but 
quadratically  integrable  functions.  The  corresponding  "interpo- 
lant"  is  simply  given  by  the  linear  combination 
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a  w 

f(x)  =  l  £m<p0(m)(x)  ■  (1.2) 

—  CO 

Note  that  the  coefficients  of  the  linear  combination  are  simply 
the  function  values  at  the  grid  points.  In  this  context  f  denotes 
the  interpolating  function,  which  is  an  approximation  to  the  orig*- 
inal  (and  unknown)  function  f,  f  =  f  ;  4>om)(x)  is  the  base 

function  corresponding  to  the  point  x  =  m  .  Observing  that 

(x  -  m)  (1.3a) 


and 


f  ■  f(m)  , 
in 

equation  (1.2)  can  alternatively  be  written  as 


(1.3b) 


f(x)  =  l  f(m)<fc0(x  -  m)  ,  (1.4) 

—  oo 


which  represents  a  product  of  the  infinite  vector  f  and  the 
vector  $o  whose  elements  are  space-shifted  base  functions. 


1.1.1  The  spectrum  of  the  step  function 


Everybody  knows  from  school  algebra  that  mathematical  oper¬ 
ations  can  be  considerably  simplified  if  performed  in  the 
"logarithmic  domain".  Multiplications  and  divisions  are  transformed 
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to  additions  and  subtractions,  etc.;  the  logarithmic  function 
In ( • )  transforms  from  the  "space  domain"  into  the  "logarith¬ 
mic  domain",  the  exponential  function  exp ( • )  transforms  in 
the  counter-direction.  Therefore,  ln(-)  ana  exp ( • )  are 
termed  "transform  pairs". 

As  far  as  functions  are  concerned,  a  particularly  useful 
tool  is  provided  by  the  Fourier  transform.  If  the  space  domain 
function  f(x)  is  absolutely  integrable, 


/  | f ( x ) ; dx  '  ®  ,  (1.5) 

—  oo 


then  the  Fourier  transform  and  its  inverse  transform  exist. 


F(r)  =  /  £(xje  i2lTr]Xdx  (1.6a) 

—  <X> 

00 

f (x)  =  /  F ( n ) e 1 4 ’  n xdn  .  (1.6b) 

-•00 

Note  that  the  Fourier  transform  F(n)  is,  in  general,  a  complex 
function.  f(x)  and  F(n)  are  frequently  called  "Fourier  trans¬ 
form  pairs".  f(x)  is  the  space  function,  F(n)  the  correspond¬ 
ing  frequency  function,  n  the  frequency. 

Let  us  consider  the  Fourier  transfrom  of  the  step-base  func¬ 
tion  <t  (x)  as  defined  by  equation  (1.1) 


J>  3  (  n  )  =  f t  (x) e  1^'rr’ xdx  = 


-  y2 


■  l  Z  rr  n  x 


dx 


Since  -?o(x)  is  symmetric,  the  transform  00(n)  is  real  and  can 
easily  be  verified  to  equal 
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This  function,  which  is  also  referred  to  as  "sine"  function  in 
literature,  is  of  particular  importance  in  connection  with 
Fourier  transforms.  It  has  zeroes  at  all  integers  n  =  ±1,  ±2,  . 
and  assumes  the  value  1  at  the  origin  n  =  0  (Fig.  1.2)  . 


In  full  analogy  to  (1.5)  and  (1.6),  the  Fourier  transform 

of  an  absolutely  summable  sequence  f  =  (f  }  , 

m 


03 


—  OO 


f 

m 


<*  oo 


is  defined  as 


(1.5)  ' 


F(n) 


L  f-e 


(1.6a)' 


a 


its  inverse  assumes  the  form  (Fuller,  1976,  p.  115  ff.) 


f 

rn 


4 

-  1/., 


F(n) 


e  i2ffrimdn 


m  =  0,  ±1 , 


(1  .6b)  ' 


In  the  following  we  show  that  the  sequence  {f  } 

IQ 

a  band-limited  process. 

According  to  equation  (1.3b)  the  sequence 
sample  of  a  continuous  function  f(x)  ,  f  =  f(m) 

IQ 

ation  functional  can  also  be  expressed  by 


represents 

{ f  }  is  a 
m 

.  This  evalu- 


00 

f  =  j  f  (x)  6  (x  -  m)  dx 

m  J 

—  oo 


with  6 ( x)  denoting  the  Dirac  5  -  distribution.  The  sequence  { f m } 
can,  consequently  be  expressed  by 


00 

( f m }  =  /  f (x) { c  (x  -  m)  }dx  , 

—  oo 

which  we  simply  write  as 


f(x)  =  fix)  •  l  i(x  -  m)  ,  (1.8) 

•oo 

where  the  product  is  to  be  interpreted  in  the  sense  of  distribu¬ 
tion  theory.  The  function 

oo 

A (x)  =  l  6  (x  -  m) 


(1.9) 
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Ls  called  an  "impulse  comb”,  its  Fourier  transform  can  be  shown 
r.o  reproduce  this  impulse  comb  (Brigham,  1974  ,  p.  19ff)  , 


A  (  n  )  = 


i  (  n  —  m) 


(1.10) 


According  to  the  convolution  theorem,  a  product  of  two  functions 
in  the  space  domain  corresponds  to  a  convolution  of  the  corre¬ 
sponding  transforms  and  vice  versa.  Therefore,  the  Fourier  trans¬ 
form  of  f (x) ,  defined  in  equation  (1.8),  is  given  by 


F  (  n  )  =  F(o)  *  A(n) 


(1.11) 


/  F (  : ) 6 ( n  -  m-  T)di  , 


which  finally  equals 


F  ( n  )  =  l  F  (  n  -  m)  . 


(1.11) 


Consequently,  F(n)  is  periodic  with  period  equal  to  1  and  the 
corresponding  frequencies  range  between  --y  <  n  <  ^  concluding 
the  end  of  the  proof  that  the  sequence  {f m }  represents  a  band- 
limited  process.  Equation  (1.6a)'  can  be  verified  by  modifying 
(1.8)  to 


?(x)  =  l  £  Six  - 
u  m 


(1.8) ’ 


Since  a  space  shift  by  m  corresponds  to  a  multiplication  by 
e  w  nm  in  the  frequency  domain,  and  since  the  transform  of 
the  Dirac  distribution  equals  1  , 


F  (p )  =  !  f  e 

*•  m 
—  00 


—  i  2  7T  nm 


(1)  ...  *  denotes  the  convolution 
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follows  immediately. 

This  result  can  be  directly  applied  to  the  problem  of 
calculating  the  Fourier  transform  of  a  step  function  as  defined 
by  equation  (1.4)  :  the  transform  of  the  space  shifted  base  func¬ 

tion  $(x  -  m)  ,  with  j(x)  defined  by  (1.1)  and  its  transform, 
by  (1.7),  equals 


simr  n 
t  n 


-  i  2  Ti  nm 

e  , 


and  consequently,  the  Fourier  transform  of  the  step  function 
(1.4)  is  given  by 


F(n) 


oo 

sinTtn  r  ,  -i2irnm 

-  )  f  e 

to  L  m 

—  OO 


(1.12) 


Let  us  split  up  F(n)  into  a  product 


F  (  n  )  =  F i ( n )  •  F; (n) 


(1.13a) 


with 


F  i  ( n ) 


sinff  n 

7T  n 


(1.13b) 


and 


F2  (n) 


-  1  2  ir  n  m 


(1.13c) 


The  following  will  be  shown  in  the  sequel: 

Every  local  interpolation  function  has  a  spectrum  of  the  form 
(:.13a)  with  Fj(n)  independent  and  ( n  >  dependent  on  the 
interpolator. 


Fjfn)  characterizes  the  interpolator  base-function  and  will, 
therefore,  be  called  the  interpolator- char acteristic;  it  acts 
as  a  frequency  dampening  function. 


1.2  THE  PIECEWISE  LINEAR  INTERPOLATION 


Linear  interpolation  is  a  generally  adopted  tool  for 
many  different  kinds  of  purposes.  As  far  as  the  representation 
of  two-dimensional  surfaces  is  concerned,  it  is  frequently 
employed  in  connection  with  a  finite  element  representation. 
Here  we  are  primarily  interested  in  the  linear  interpolation 
of  an  infinite  sequence  of  data  located  at  the  cardinal  numbers 
and  in  its  spectrum. 

The  base  function  is  the  continuous  and  strictly  local 
function  (Fig.  1.3) 

ll  -  ( x |  for  |x|  s  1 

« l (x)  =  |  (1.14) 

jo  else 


It  is  an  element  out  of  K1  ,  the  space  of  continuous  functions 
having  quadratically  integrable  first  order  derivatives.  Its 
support  has  length  2  ,  this  is  1  unit  more  than  that  of  the 
step  base-function.  The  corresponding  interpolant  is  simply 
given  by  the  linear  combination 


=  l  f  <C  i  (x  -  m)  . 
^  m  ‘ 


f  (x) 


(1.15) 
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1.2.1  The  spectrum  of  the  piecewise  linear  interpolant 


The  Fourier  transform  of  the  base  function  (1.14)  is 
given  by 

-»!  (n)  =  /  (1  -  ;xi)e'i2,nxdx 

-  l 

which  can  easily  be  shown  to  equal 

,  ,  ,  _  fsimm 

l  *n 

Since  i>  -  is  the  result  of  the  convolution  between  $>0  and 
*9  ,  it  follows  from  the  convolution  theorem  that  the  correspon¬ 

ding  Fourier  transform  is  a  simple  product 


(1.17) 


t ;  —  $9  •  $9  t 


which  is  naturally  fulfilled  by  equation  (1.17)  with  4>0  given  by 
(1.7)  . 

Taking  into  account  (1.13a,c)  and  (1.17),  the  Fourier  trans¬ 
form  of  the  piecewise  linear  and  continuous  interpolant  is  simply 
given  by 


F(n)  =  !'5iSHLU  l  f  e“  i  2  t  n  m  _  (1.l8) 

-  n  1  L  m 

s  '  —  TO 

The  piecewise  linear  and  continuous  interpolation  function 
can  be  differentiated  along  the  real  line  excluding  the  cardinal 
numbers : 


a 

D  f  (x)  =  7  f  D  $  i  (x  -  m) 

x  '  *  m  x 

—  on 


with 


Dx$ : (x)  = 


for  -1  <  x  <  0 


0  <  x  <  1 


(1.19) 


else 


In  order  to  find  the  spectrum  of  Dxf  we  proceed  as 
follows:  the  interpolant  f  equals  the  inverse  Fourier  trans¬ 
form 


f  ( x )  =  /  F(  n)  el2TTnxdn  ; 


we  differentiate  f  and  obtain 


D  f  ( x )  =  i2rr  /  F(n)el27Tnxndn  ; 

v  4 


denoting  the  spectrum  of  Df  by  G  ,  we  find 


G ( n )  =  /  D  f(x)e 


-  i  2  n  n  x 


.  j  -  i  2  it  i  x  . 


=  i2:i  /  [  /  F(n)ei  TTrlxndnie 


We  integrate  first  with  respect  ot  x  ; 


/  ei2’Tx(""lJ)dx  =  6  ( n  -  u) 


according  to  the  theory  of  distributions  (Brigham,  1974,  p.  230) 
consequently , 


G(n)  =  12t7  /  F  ( n )  n  <5  ( n  - 


u )  dn 
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and  considering  the  properties  of  the  Dirac  distribution,  we 
obtain  for  the  spectrum  of  the  differentiated  piecewise  linear 
and  continuous  function 


G  ( -i )  =  i  2  t  r  F  (  n  ) 


(1.20) 


1  .3  THE  QUADRATIC  SPLINE  INTERPOLATION 

Both  the  base  function  $0  for  the  step  function  "inter¬ 
polation"  and  $ i  for  the  piecewise  linear  interpolation  are 
such  that  its  support  does  not  exceed  two  units.  This  is  the  rea¬ 
son  why  the  coefficients  of  the  linear  combination  are  identical 
with  the  function  values  at  the  grid  points.  This  is  no  longer 
rrue  for  a  quadratic  spline  base  function  which  has  a  support 
of  3  units.  At  the  end  of  section  1.2  we  have  briefly  stated 
that  $  i  is  the  result  of  a  convolution  between  $9  and  . 

It  is  left  as  an  exercise  for  the  reader  that 

(j  2  =  4>  i  *  <t>  g  (1.21) 


with 


$2  (X) 


j  l  (-1)  j(,3)  (X  +  |  -  j  )2 

*  3-0  J 


(1.22) 


where 


for  x  >  0 
else 

Explicitly  written  »2 (x)  is  as  follows 

1_ 

2 

x|  S  |  (1.22)  ' 


It  is  evident  that  the  support  of  <p2  is  the  open  interval 
3  3 

(--  ,  j)  .  Of  course,  $2  is  symmetric  and  assumes  a  maximum 
value  of  j  at  the  origin,  $2(0)  =  j  .  It  can  be  easily  verified 
that  ^ 2 ( x )  is  continuous  and  in  addition,  continuously  differen¬ 
tiable.  (Note  that  its  first  derivatives  vanish  at  x  =  -~  and 
x  =  j  . )  It  is,  therefore,  an  element  of  K2  ,  the  space  of  con¬ 
tinuous  and  continuously  differentiable  functions  having  quadrat- 
ically  integrable  second  order  derivatives.  Its  integral  is  one, 


2>  2  ( x )  = 


f  2,3 

i  "X  4 

1  fl  _  |x,'. 

2  2  1  1  I 


for 


! x  i 


for  j  £ 


else 


(x) 


1  . 


Obviously  <J>2  is  not  a  sampling  function  (i2(k)  /  5  ) 

which  can  be  directly  seen  from  Fig.  1.4.  Therefore,  it  is  not 
possible  to  express  the  interpolation  function  in  a  form  like 
(1.4)  or  (1.15).  With  other  words,  the  coefficients  of  the  linear 
combination  (g  } 


FIG.  1.4  Base  function  for 


f(x)  =  I  gm$2 (x  -  m)  (1.23) 

—  oo 


do  no  longer  equal  the  function  values.  They  can,  however,  be 
pressed  as  a  linear  combination  of  the  function  values.  This 
shown  as  follows. 

Let  us  introduce  a  cardinal  spline  or  sampling  spline 
02(x)  such  that 


la 


9 


Since  the  elements  of  A  ,  a.^  ,  depend  only  on  the  difference 
k-j  ,  it  follows  that  A  is  a  Toeplitz  matrix  of  infinite 
dimension  and,  therefore,  equivalent  to  a  circulant  matrix.  Such 
kinds  of  systems  are  particularly  easy  to  solve  employing  Fourier 
transform  methods  (Grenander  and  Szego,  1958)  :  If  a ( A }  and 

e(')  denote  the  discrete  Fourier  transforms  of  any  row  (or 
column)  of  A  and  the  vector  e  ,  respectively,  the  solution 
vector  u  is  simply  given  by  the  inverse  Fourier  transform  of 
e(»)  /  a(\)  .  (This  is  basically  due  to  the  fact  that  (1.27) 

represents  a  discrete  convolution  v/hich  reduces  to  a  correspond¬ 
ing  product  in  the  frequency  domain.) 

a  (A)  =  rr  ( e  i'  +  6  +  e1')  =  j(3  +  cos  A)  , 

5(a)  =  i  • 

The  inverse  Fourier  transform  of  e(A)  /  a ( x )  is  simply 


a.  =  — 

K 


(  coskA 
j  3+cos'' 


dA 


which  can  be  shown  to  equal  (Gradshteyn,  1971,  p.  380,  No.  3.613) 


a,  =  /2  (-3  +  2/2)  ^ki  .  (1.29) 

k  -  k 

With  (1.29)  the  cardinal  spline  l-2  ,  shown  in  Fig.  1.5,  can  be 
expressed  by 


'll  2 


>1  l  (-3  +  2/2) 


(x  -  k) 


(1.26)  • 


This  infinite  sum  can  be  reduced  to  a  finite  sum  consisting  of 


only  3  terms  if  we  take  into  account  the  properties  of  the 

base  function  $2  ;  $2  ^as  a  compact  support  which  is  the  open 

3  3 

interval  (~j  ,  j)  (see  Fig.  1.4),  $2  is  symmetric  (?>2(-x)  = 

=  i i2(x)).  Observing  the  symmetry  of  a.  ,  a  =  a  ,  equation 

K  “*  K  K. 

(1.26)'  reduces  to 

k ,+  2 

C-  ;  ( X )  =  l  « k  2  (  i  x  |  “  k ) 

k  =k  1 

with  (1.26)" 

k j  =  int ( |  x I  +  0.5)  -  1 

The  verification  that  v2 (x)  is  in  fact  a  sampling  function 
fulfilling  equation  (1.24),  is  left  as  an  exercise  for  the  reader 


FIG.  1.5  Cardinal  spline  >,(x) 
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It  seems  to  be  worthwile  to  point  out  that  ■.  has  no  Icnqer  a 
finite  support.  This  is  somewhat  astonishing  since  the  cardinal 
sampling  function  of  the  space  K1  ,  denoted  by  p  ,  had  still 
a  finite  support.  This  is  basically  due  to  the  fact  that  the 
inverse  of  the  tri-d.iagonal  transformation  matrix  A  is  a  full 
infinite  matrix  for  K2  ;  the  corresponding  matrices  for  K: 
and  K-  are  identity  matrices  :  the  base  functions  and 

fulfil  equation  (1.24)  and,  therefore,  these  base  functions 
are  already  sampling  functions.  Furthermore,  the  sum  of  all 
'  s  equals  1 


7  1,  =  /2  : 1  +  2  7  (-3  +  2/2) k]  =  1  ,  (1.30) 

-«  l 

and,  therefore,  considering  the  fact  that  the  integral  of  i? 
is  unity,  also  the  integral  of  >i»2  equals  1  , 

co 

/  'J/2  (x)  dx  =  1  . 


1.3.1  The  spectrum  of  the  quadratic  spline  interpolant 


The  Fourier  transform  of  the  quadratic  spline  interpolant 
( equation  (1.25))  , 

A  °° 

f  (x)  =  7  fn-ji?  (x  -  m) 
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is  given  by 


F(n)=  If  J  v  2  ( x  -  m) 


-  i  2  n  n  x  . 


which,  according  to  the  translation  property  of  Fourier  trans¬ 
forms  reduces  to 


F(n)  =  /  v2(x)e“i2Trnxdx  •  \  f  e‘i27Tnm 


Considering  the  definition  of  the  cardinal  sampling  function 
given  by  equation  (1.26}  and  taking  into  account  the  symmetry 
of  v 2  >  the  above  integral  can  be  written  as 


/  '(j2(x)e  1",rnxdx  =  St  ( n )  i  a  cos  ( 2-rnk) 


1.31) 


Here  •}  2  ( t )  denotes  the  Fourier  transform  of  the  base  function 
(x)  .  Since  is  the  result  of  a  convolution  between  <f ; 

and  (equation  1.21)  ,  it  follows  from  the  convolution  theorem 

that  i ^  is  a  product  of  and  ;  <t>  i  in  turn  is  a  product 

of  and  with  given  by  equation  (1.7)  .  Therefore, 


P  2  (  n )  =  P  n  ) 


(1.32) 


The  infinite  sum  in  equation  (1.31)  needs  still  to  be  calculated. 
Since  a  and  the  cosine  function  are  symmetric  and  moreover, 
according  to  equation  (1.29) 


=  »  2 


=  -3  +  2,2 
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we  can  write 


D  JO 

»kcos2-Tnk  =.'2(1  +  2  T  akcos2mk)  , 
—  00  1 


and  due  to  (Gradshteyn,  1971,  p.54,  No.  1.447)  the  infinite  sum 
can  be  expressed  in  a  closed  form  such  as 


oo  1  —  ^ 

y  a,  cos2-r-,k  =  /2  -  .  (1.33) 

^  K. 

-«  l-2aCOs2ir  n+a  ^ 

With  a  as  defined  above  this  equation  can,  after  trivial  algeb¬ 
raic  operations,  be  reduced  to 


—  oo 


a  cos  2  it  nk 

K 


4 

3  +  cos2tti 


(1.33)1 


With  (1.31)  we  finally  obtain  the  Fourier  transform  of  the 
quadratic  spline  interpolant 


F(n) 


4 

3+cos2tt  n 


(  \  •>  » 

I  sinTr n  J  y  f  -i2Trnra 

i  ^  ®e 


(1.34) 


which  is  again  of  the  form  (l.l 3a). 

The  quadratic  spline  is  continuous  and  continuously  differ¬ 
entiable  with  quadratically  integrable  second  order  derivatives. 
The  spectrum  of  its  first  and  second  order  derivative  can  be 
derived  in  complete  analogy  to  subsection  1.2.1.  We  obtain  for 

A 

the  spectrum  of  D*£(x) 


G i  ( n )  =  i  2  tt  n  F  ( n ) 


and  for  the  spectrum  of  D^f(x) 


(1 .35a) 
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It  is  obvious  that  the  support  of  c  j  is  the  open  interval 

{-2,2)  ,  its  length  is,  therefore,  equal  to  4  .  Of  course, 

2 

5  is  symmetric  and  assumes  a  maximum  value  of  ■=•  at  the 
origin,  i (0)  =  —  .  It  can  be  easily  verified  that  $  ,  is 
continuous  and  twice  continuously  differentiable.  (Note  that 
r,  ,  its  first  and  second  order  derivatives  vanish  at  x  =  -2 
and  x  =  2  .)  Therefore,  $3  is  an  element  of  K3  ,  the  space 
of  continuous  and  twice  continuously  differentiable  functions 
having  quadraticaily  integrable  third  order  derivatives.  Its 
integral  is  one,  03  is  a  base  function  for  the  cubic  spline 
interpolation,  but  unfortunately  not  a  sampling  function 
(i3(k)  d  3  )  which  can  be  seen  from  Fig.  1.6.  Therefore,  it 

is  not  possible  to  express  the  interpolation  function  in  a  form 
like  (1.4)  or  (1.15);  we  are  in  a  similar  situation  as  we  have 
been  with  t2  1  the  coefficents  of  the  linear  combination  are 
no  longer  the  function  values  themselves;  however,  they  can  be 
expressed  as  linear  combinations  of  the  function  values  at  the 
cardinal  numbers. 

The  derivation  of  the  cardinal  spline  or  sampling  spline 
:•  3  is  practically  identical  to  the  derivation  of  •!/ 2  as  presen¬ 
ted  in  section  1.3.  Therefore,  we  point  out  only  the  important 
steps . 

fulfilling  equation 

(1.38' 

The  infinite  vector  of  coefficients  j  1  can  be  determined  from 
the  known  properties  of  the  base  function  >3  and  the  sampling 
properties  of  : 1  postulated  by 

v,  fk)  =  ■*  .  (1.39) 

k  ,  0 


The  cardinal  cubic  spline  y3(x)  , 

(1.39)  ,  is  defined  by 

v -  ( x )  -  r  ; i (x  -  k) 


Evaluating  (1.38)  at  all  cardinal  numbers  and  observing  (1.39)  , 
we  obtain  a  linear  system  of  infinite  dimension. 


A«  =  e  , 


with  i  and  e  denoting  the  infinite  vector  of  coefficients 
{  i ,  }  and  e  the  infinite  vector  given  by  equation  (1.28). 

0 3  has  finite  support  of  length  equal  to  4  (see  Fig.  1.6)  and 
does  not  vanish  at  the  cardinal  numbers  x  =  -1,  0,  1  .  There¬ 
fore,  the  symmetric  matrix  A  has  only  3  non-vanishing  dia¬ 
gonals.  The  main  diagonal  elements  are  constant  and  equal  to 
o5(0)  =  —  ,  the  elements  of  the  other  two  diagonals  are  also 

constant  and  equal  to  « 3  ( —  1 )  =  4>3(1)  =  i  , 


(1.40) 

The  Elements  of  A  ,  a  .  ,  depend  only  on  the  difference  k  -  j  ; 

]  k 

therefore,  A  is  a  Toeplitz  matrix  of  infinite  dimension  and 
equivalent  to  a  circulant  matrix.  The  solution  vector  is  found  by 
Fourier  transform  methods  in  the  very  same  way  as  it  was  obtained 
in  section  1.3  :  the  Fourier  transform  of  a  row  (or  column)  of 
A  is  given  by 


a(A)  =  i  +  4  +  e1*]  =  J  ( 2  +  cosM 

The  solution  vector  a  is  obtained  by  an  inverse  Fourier  trans- 


28 


form  of  e(\)/a(>)  with  e(.\)  =  1  as  in  the  foregoing  section, 


a 


k 


3 


cos  k  • 
2+cos  > 


d> 


which  can  be  shown  (Gradshteyn,  1971,  p.  380,  No.  3.613)  to 
equal 

_  _  1  k  ' 

a  =  a  =  /3  (-2  +  /3) '  '  .  (1.41) 

k  -k 

With  (1.41)  the  cardinal  cubic  spline  ^  3  ,  shown  in  Fig.  1.7, 
can  be  expressed  by 

co 

*3(x)  =  /3  l  (-2  +  >3)  ' *  1 >  3 (x  -  k)  .  (1.38)’ 

—  00 

The  infinite  sum  reduces  to  a  finite  one  of  only  4  elements, 
when  taking  into  account  the  properties  of  the  base  function  .  3  : 


kj  +  3 

4-)(x)  =  T.»«3(!x!-k)  (1.38)" 

k  =  k  i 


with 


k ;  =  i nt ( j  x 1 )  -  1 

The  reader  is  invited  to  verify  that  v  ->  ( x )  is  in  fact  a  sampling 
function  fulfilling  equation  (1.39).  The  reasoning  for  the  unlimi¬ 
ted  support  of  >1 3  can  be  carried  over  from  section  1.3  without 
any  change. 

This  sum  of  all  coefficients  can  be  verified  to  equal  1, 


FIG.  1.7  Cardinal  spline  1 3 ( x ) 


y  a  =  /3  [1  +  2  l  (-2  +  /3)k]  =  1  ;  (1.42) 

-»  *  l 

since  the  integral  of  >3  is  unity,  it  follows  from  the  definition 
of  >  (cf.  equation  (1.38))  that  the  integral  of  p  3  is  also  1 

I  -j  3  (x)  dx  =  1 
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1.4.1  The  spectrum  of  the  cubic  spline  interpolant 


The  Fourier  transform  of  the  cubic  spline  interpolant 


f  ( x )  =  ■  f  v  i  ( x  -  m ) 

u  m 


is  given  by 


f  /  ^3(x  -  m)e  l2TTnxdx  , 


m  ' 

“CO  —  CO 


which,  according  to  the  translation  property  of  Fourier  transforms, 
reduces  to 


F(n)  =  '4j  3  (x)  e~i2lTnxdx  •  \  f  e"i27rnx, 

L  m 

—  03  —  TO 


Considering  the  definition  of  the  cardinal  sampling  function  v  3 
given  by  equation  (1.38)  and  taking  into  account  the  symmetry  of 
v-  ,  the  above  integral  can  be  written  as 


/  *3  (x)e  li7Tnxdx  =  4  5  ( n )  \  1  cos(2wnk) 


(1.43) 


Here  t>  3  (  n )  denotes  the  Fourier  transform  of  the  base  function 
: - (x)  .  Since  33(x)  is  the  result  of  a  convolution  between 

and  (equation  1.36)  ,  it  follows  from  the  convolution  theorem 

that  j3  is  a  product  of  and  i0  ;  with  (1.7)  and  (1.32) 

1  is  simply  given  by 


With  (1.43)  we  finally  obtain  the  Fourier  transform  of  the  cubic 
spline  interpolant 


32 


F(n) 


sinn  n 


2  +  COS2nr,  nn  J  - 


y  f  e 

m 


■12nnm 


which  is  again  of  the  form  (1.13a). 


(1.46) 


The  cubic  spline  is  continuous  and  twice  continuously 
differentiable  with  quadratically  integrable  third  order  deriv¬ 
atives.  The  spectra  of  these  spline  derivatives  can  be  obtained 
in  complete  analogy  to  subsection  1.2.1.  We  obtain  for  the 
spectrum  of  D^f(x) 

Gk(n)  =  (i2nn)kF(n)  ,  k=0,  ...,  3  .  (1.47) 


1.5  HIGHER  AND  HIGHEST  ORDER  INTERPOLANTS 


In  the  foregoing  sections  we  have  been  dealing  with  low 
order  interpolants  on  the  sequence  of  cardinal  numbers.  All  these 
functions  turned  out  to  be  closely  related  to  each  other.  Before 
we  present  a  condensed  form  of  the  whole  array  of  relations,  we 
would  like  to  draw  attention  to  interpolation  functions  of  higher 
order  continuous  differentiability;  we  shall  end  up  with  the 
highest  order  interpolant  possible. 

It  was  shown  in  the  last  sections  that  the  base  function 
for  the  space  KJ  is  the  step  function  -pQ  as  defined  by  (1.1), 
and  that  the  base  functions  for  higher  order  spaces  result  from 
convolutions  of  pQ's  .  In  general,  the  base  function  $  for 
the  space  Kn  is  given  by 


33 


combination  of  base  functions  <j>  , 

n 


v  (x)  =  7  a  ®  (x  -  k) 

fl  —CO  *  r* 


(i .53; 


The  infinite  vector  >a  }  is  determined  from  the  sampling  property 

k 

of  V 

n 

•;>  ( k )  =  .5  ,  (1.54) 

n  k  ,  0 

and  the  properties  of  3  ;  it  leads  to  the  solution  of  the  in- 

n 

finite  linear  system 


A,i  =  e 


(1.55) 


with 


and 


■  •  •  x  -i  ,  rx  0  ,  a  1  ,  .  .  .  , 

.*•0  ,  1  /  Qf  ...'  . 


The  matrix  A  has  very  special  properties:  it  is  an  infinite 
symmetric  matrix  of  band  structure;  all  diagonals  have  identical 
elements.  Such  a  matrix  is  known  as  a  Toeplitz  matrix  which  is 
equivalent  to  a  circulant  matrix  due  to  its  infinite  dimension. 

The  bandwidth  N  (  =  number  of  non-vanishing  diagonals) 
depends  only  on  the  support  of  the  base  function  j>  : 


if 

n 

is 

even 

if 

n 

is 

odd 

E.g.  :  For  the  linear  interpolation  (n  =  1)  ,  A  is  a  pure  dia¬ 

gonal  matrix,  moreover,  it  is  the  unit  matrix,  therefore 
and,  as  a  consequence  (equation  (1.53))  ,  the  original  base 


e 
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function  i  •  and  the  corresponding  sampling  function 


cide . 


Denoting  a  row  (or  column)  of  A  by  a  ,  a  =  ;  a,  '•  ,  the 


elements  { a,  }  are  given  by 


a,  =  t  (k) 
k  r. 


(1.57) 


The  inverse  of  A  can  be  found  via  Fourier  transform  methods 
Denoting  the  discrete  Fourier  transform  of  a  by  a  , 


)  (k)cosk\ 
L  n 


int (n/2) 

$  (0)  +  2  7  p  (k)coskA  , 

n  ^  n 


(1.58: 


and  by 


r  -  i  k  A  , 

L  e  ue  =  1  / 


the  coefficient  vector  of  unknowns  (a,  1  can  be  directly  founc 

k  u 


by  the  inverse  Fourier  tranform 


"k  2- 


e  (  a  )  _  i  k  •  j  . 

— -  e  d  \ 


i  r  i 


cos  k'd  i 


(1.59) 


Chis  integral  plays  a  central  role  within  the  framework  of  splines 


if  only  1  frequency  is  represented  in  a  ( '•  )  (apart  from  the 
zero  frequency)  ,  equation  (1.59)  admits  a  closed  expression 
and  the  a^s  can  be  found  directly  (n  -  3)  ;  if  more  than  i 
frequency  is  represented  (n  •  3)  ,  this  simplicity  is  no  longer 

present . 

The  cardinal  spline  of  degree  n  ,  „  ,  as  defined  by 

n 

equation  (1.53)  ,  has  unlimited  support  for  n  >  1  ;  the  higher 
n  ,  the  weaker  is  its  tendency  to  tend  to  zero  with  increasing 
argument,  and  therefore,  the  less  is  its  local  behavior.  Its 
spectrum  has  the  general  form 


F 

n 


(  n ) 


H  ( n ) 
n 


s in^  n • 

,  tt  r,  j 


CO 


l 

■—  oo 


f  e 
m 


2  i  n  m 


(1.60) 


where  H  (  r. ) 

n 

H__  1  ,  for 

Explicitly 


is  characteristic  for  the  order  n  .  For  n  =  0,  1 
n  =  2,  3  H  n  is  given  by  (1.33)'  and  (1.45). 


H  (-,) 
n 


)  t  COS2irr|k 
-  k 

—  » 


(1.61) 


with  i  ;  denoting  the  vector  of  coefficients  for  the  degree 

k 

n  in  consideration.  These  functions  H  (c)  show  the  very  inter¬ 
im 

estmg  property  of  compensating  the  dampening  function 

r*  +  I 

(sinm/rn)  ‘  better  and  better  with  increasing  n  . 

Let  us  put  now  the  obvious  question: 

:v.*,  a  c  kind  of  function  is  the  limit  csrihnal  spline 


limy  ( x  ) 
n 
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Obvious iy 


f  (x)  =  "  f  v  (x  -  n) 

»  m  * 

should  interpolate  the  data  vector  (f  }  which  is  supposedly 
an  element  of  1*  ,  the  space  of  quadratical ly  summable  finite 
differences  of  all  orders  k  , 


k=o 

The  cardinal  function  itself  should  correspondingly  be  an  element 
of  K°°  ,  the  space  of  quadratically  integrable  derivatives  of 
all  orders  k  , 


00  k 

U  K*  . 

k  =0 


There  is  a  one-to-one  correspondence  between  the  spaces 
K?  (Schoenberg,  1973).  The  unit-sequence  (...,  0,  i,  0, 
surely  an  element  of  12  .  It  can  be  shown  (Schoenberg, 

that  the  function 


* 


. . . )  is 
1973) 


sinirx 

TX 


is  the  only  element  of  K2  which  interpolates  the  unit  sequence. 
Therefore,  we  conclude  that 

, .  ,  ,  sinrx 

lim<j/  (x)  =  - 

n  ttx 

n— ” 

The  corresponding  limit  interpolation  function  is 


(1.62) 


can  be  said  to  be  the  cardinal  spline  of  highest  possible 
degree  (»)  ,  f  the  corresponding  interpolant.  Its  tendency 

to  approach  zero  is  rather  weak  and,  therefore,  its  interpolat 
behavior  cannot  be  considered  local .  (Fig.  1.8). 


FIG.  1.8  Cardinal  spline  .^(x) 

The  spectrum  of  •*»  can  easily  be  shown  to  be  a  unit 

on 

step  function.  Actually  it  has  already  been  shown  in  subsection 
1.1.1  that  the  spectrum  of  the  unit  step  function  equals  sin 
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T  (n) 

a 


OO 


\  sirvrx 

J  -x 


’urnx 


dx 


1  for 
y  for 
0  else 


1_ 

2 

1 

2 


(1.64) 


(3raashteyn,  1971,  p.  428,  No.  3.74 
all  n  '  -*  •=■  ,  we  conclude  that 

Z  +■ 

interpolation  function  fs  does  not 
than  "  >  i  .  (1.64)  and  the 

into  account ,  the  Fourier  transform 


1) .  Since  V  (n)  =  0  for 

on 

is  a  low  pass  filter  -  the 
contain  frequencies  higher 
translation  property  taken 
of  f  is  given  by 


F  (n) 


-  V  f  e'i2rrun  for 

m 

( 

1  0  for 


(1.65) 


(Here  we  did  not  consider  the  case  in1  =  i  .)  The  relation  between 
and  V ,  is  quite  remarkable  ar.d  deserves  special  attention: 


The  cardinal  spline  of  highest  possible  degree  equals  the 
spectrum  of  the  cardinal  spline  of  lowest  possible  degree  and  vice 
versa  . 


=  f 


(1.66) 
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1.5.1  Guidelines  for  cardinal  spline  interpolation  of  arbitrary 
degree 


The  detailed  and  straightforward  presentation  in  sections 
1.1  to  1.4  will  be  complemented  by  some  additional  and  quite  use¬ 
ful  relations;  the  procedures  involved  starting  from  the  choice 
of  an  appropriate  base  function  and  ending  w_th  the  interpolation 
formula  are  presented  in  a  flowchart-like  diagram. 


a)  The  first  additional  information  is  related  to  the  function 

values  of  a  base  function  ;  at  the  cardinal  numbers.  Remember 

n 

that  they  represent  the  entries  of  the  transformation  matrix  A 

(equations  (1.55)  and  (1.57))  which,  in  turn,  provides  the  link 

between  the  base  function  and  the  corresponding  sampling  function, 

therefore,  they  simply  need  to  be  known. 

Basically  there  are  two  approaches,  a  direct  one  which  derives 

the  function  values  employing  Peano's  theorem  (Schoenberg,  1973, 

p.  2  ff.)  and  an  indirect  one,  which,  according  to  the  author, 

is  simpler  and  more  transparent: 

We  know  from  (1.51)  that  the  Fourier  transform  of  the  base 

function  t>  is  the  (n+l)th  power  of  sinin/^n  , 
n 


r  .  (  n ) 


!  \ntl 

!  5  int!  r\  \ 

i  t  n  I 


The  base  function  t  ,  therefore,  is  obtained  by  the  inverse 

n 

Fourier  transform  of  $  , 

n 


;  (x)  =  2  1  't  ( n )  cos2 t  nxdn 

n  J  n 

0 


Consequently,  the  values  of  the  base  function  at  the  cardinal 
numbers  is  simply  obtained  by  replacing  the  variable  x  by 
an  integer  k  ( 

V) 

t  (k)  =  2  ■  — n”1 j ‘  cos2"kndn  .  (1.67) 

n  <  i  *  n  J 

3 

With  (Gradshteyn,  1971,  p.  471,  No.  3.836)  this  integral  can 
be  expressed  as 


int  A  <-  1  k  ' 


_1_ 

n! 


5=0 


(k)  = 


(-11 J '"j1] * 


k 


|  r1 

-  j  ‘  for  0  ;  k  ,  in 


for  ’ k '  •  int 


n  i 


12, 


(1.68) 


It  can  easily  be  shown  to  be  equivalent  to  an  expression 
obtained  by  Schoenberg  (1973) 


n+  l 

>  ( k)  =  ~r  l  (-1) 


1  “r  * /  , , j ' n+ll  | n+1 


n : 


j  =0 


■  +  k  -  j  I 

D  )  [  2  J  i 


inking  into  account  the  definition  of 


x  '  for  x  >  0 


x  = 


(1.69) 


0 


else 
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Either  (1.68)  or  (1.69)  allows  us  to  verify  the  identity 


k 


oo 

l 


i  (k) 
n 


=  1 


V  n  €  N  q 


(1.70) 


b)  The  second  additional  information  concerns  the  calculation  of 
the  cardinal  spline  sampling  functions  from  the  corresponding 
base  functions.  In  all  foregoing  sections  we  have  presented 
one  way  of  obtaining  the  sampling  functions  which  can  be 
sketched  as  follows: 


{ a  }  =  { ®  (k)  } 
k  n 

i 

oo 

a  ( n)  =  y  a,  cos2irko 

k 

—  OO 

f  T 

i  =  2  ^  rr— —  cos2nkndn 

k  ,  J  a('l} 

Vx)  =  l  Vn(x  ~  k)  • 

_  oo 

The  other  way  of  obtaining  'J>  avoids  the  explicit 
calculation  of  the  infinite  vector  a  =  {a.  }  ;  this  and  the 
above  (at  least  theoretically)  infinite  summation  can  be  repla¬ 
ced  by  an  inverse  Fourier  transform.  Observe  that  the  Fourier 

transform  of  ■„  is  given  bv 

n 


v  (  n ) 

n 


o 

i  (  n )  /  cos  2  ik",  ; 

n  k 


according  to  (1.55) 


which  transforms  to 
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30  CD  >  _  , 

a,  cos2’Tkn  =  !  >  a,  cos2i!;n  :  =  a- !  (  n  )  : 

-  k  lb*  l 

(note  that  the  vector  a  represents  a  row  or  column  of  A" 1 ) . 
Taking  this  relation  into  account,  the  Fourier  transform  of 
V  is  given  by 


y  „  ( i ) 


*n  <  n) 
a  ( i ) 


(1.71) 


and  the  cardinal  spline  of  degree  n  by  its  inverse  Fourier 
transform 


<l>  (x) 
n 


r  n  + 

isinm 

l  m  ) 


cos2n  nxdn 


'  4>  (k)  cos2’rkn 
n 

—  00 


(1.71)  ' 


We  have  thus  avoided  the  actual  calculation  of  the  coefficients 
(i  }  .  This  representation  shows  very  clearly  that  only  the 
basis  spline's  function  values  at  the  cardinal  numbers  are  re¬ 
quired  to  obtain  the  corresponding  cardinal  spline;  it  has, 
however,  no  practical  consequences , since  it  is  much  easier  to 
perform  a  Fourier  transform  (for  the  cardinal  numbers  only) 
once  for  all,  and  then  to  combine  the  resulting  ,1k's  linearly 
with  the  base  functions  as  in  (1.53) .  We  will  need  (1.71)  for 
a  different  reason  in  the  next  chapter. 

The  following  graphs  illustrate  quite  clearly  the  relations 
between  spline  base  functions,  their  Fourier  transforms,  and  cor¬ 
responding  cardinal  splines  and  their  Fourier  transforms. 


9  Spline  base  functions 


of  degree  n 


10  Fourier  transforms  $ 

n 

of  degree  n 


of  spline  base  functions 


SPLINES  AND  THE  GAUSSIAN  FUNCTION 


2  . 


In  the  last  chapter  we  intended  to  introduce  the  family 
of  cardinal  splines  in  a  stepwise  direct  approach.  We  saw  that 
tnere  is  a  whole  bunch  of  relations  between  these  interpolation 
functions;  section  1.5.  focused  on  the  goal  of  finding  the  cardi¬ 
nal  spline  of  highest  possible  degree;  it  turned  out  to  be  the 
familiar  simx/irx  function. 

In  this  chapter  it  should  be  investigated  if  there  is  any 
relation  between  splines  and  the  Gaussian  bell-shaped  function. 

The  basic  starting  point  is  again  the  cardinal  spline  base  func¬ 
tion  of  lowest  possible  degree,  the  zero  degree  unit-step  func¬ 
tion  4>0  .  We  noted  already  during  the  discussion  of  higher  than 
zero  order  spline  base  functions  a  couple  of  significant  features: 

a)  the  base  function  <j> 0  has  a  constant  height  equal  to  1  and 
a  support  length  of  1  ; 

bj  the  base  function  $  can  be  expressed  as  an  n-fold  convolu¬ 
tion  of  the  unit-step  zero  degree  cardinal  spline  base  func¬ 
tion  expressed  by  equation  (1.48); 

c)  with  each  increase  by  one  degree  the  continuous  differentiability 

increases  by  one  order;  0  is  (n-1)  -  times  continuously 

n 

differentiable; 

d)  with  each  increase  by  one  degree  the  support  length  increased 

by  one  unit;  has  a  support  of  (n+1)  units;  at  the  same 

time  the  maximum  value  0  (0)  becomes  smaller;  the  functions 

n 

flatten  out  with  increasing  degree; 


I 
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c'  the  correlation  length  increases  with  increasing  degree  as 
can  be  seen  from  Fig.  1.9  ,- 

f)  all  base  functions  $  have  a  constant  integral  area  equal 
to  1  ; 

a)  the  sun  of  all  function  values  at  the  cardinal  numbers  equals 
1  for  any  r 

n 

We  are  making  now  a  side-step  and  recall  some  facts  known 
from  probability  theory:  According  to  the  central  limit  theorem, 

the  density  function  of  a  sum  of  uniformly  distributed  independent 
random  variables  tends  to  a  Gaussian  function  (Papoulis,  1965, 
p.  267),  More  specifically:  if  n  independent  random  variables 
x<  have  a  constant  density  function  which  equals  the  step  function 
:i  (equation  (1.1)),  it  follows  immediately  that 

S' x}  =  E{  >  X  :  ■  E  x  }  =  0 

1  7  1 

1  1 

The  variance  of  a  single  random  variable  is  accordingly  given  by 


E*  x-  -  =  x2dx  =  -t4r 

i  12 

-  2 

Sinc°  the  random  variables  are  independent  as  postulated  above, 
the  variance  of  the  sum  equals  the  sum  of  the  variances, 

E-  (  7  x  ) ; }  =  y  E{xz  }  ; 

r  i  .  i 

i  i 

it  follows  that  for  the  case  discussed  here  (i  =  l,  . . . ,  n) 
the  variance  is  given  by 
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n 

a2  =  l  E{x2}  =  i  (2.1) 

i  =  1 

The  density  function  of  the  sum  equals  the  (n-l)fold  convolution 
of  the  individual  densities  (which  were  assumed  to  be  equal). 


Y  (x)  =  $q(x)  *  <(>  0  (x)  *  ...  *  $0  (x) 


(2.2) 


(n-1)  times 

The  corresponding  Gaussian  density  function  is 


Y  (X) 


- -  e 

6  /2n 


and  with  (2.1) 


Y  (x) 


(2.3) 


With  7  and  y 

n  n 

approaches  yr  if 


the  central  limit  theorem  says  that 
n  increases  , 


Y 


n 


Y 


n 


Y 

n 


(2.4a) 


and  under  very  general  conditions, 

1  im  y  =  Y 
r  n 

n  -*■» 


(2.4b) 


49 


Let  us  now  come  back  to  the  spline  base  functions.  Comparing 
,1.48)  with  (2.2)  ,  we  conclude  that  the  cardinal  spline  base  func¬ 
tion  of  degree  n  ,  ■i  ,  eouals  the  density  function  of 

n  r.+  \ 

n+l  independent  random  variables  having  equal  density  functions 


5 

n 


-  Y 


r.+  1 


(2.5) 


According  to  (2.4a)  and  (2.4b),  the  Gaussian  function  C 

n+l 

defined  in  (2.3) ,  can  be  considered  an  approximation  to  * 

n 


Y 


n+  1 


0 


r. 


and  in  the  limit 


lim 

n-+°° 


V 

n+  1 


* 

00 


(2.6a) 


(2.6b) 


Therefore,  the  cardinal  spline  base  function  of  degree  n  can 
be  approximated  by  a  Gaussian  function. 


;  ( x) 


_  6*1 
/  6  n+l 

'/  it  ( n+ 1 ) 


(2.7) 


The  following  graphs  provide  a  comparison  between  exact  cardinal 
spline  base  functions  and  their  corresponding  approximations  by 
a  Gaussian  curve. 
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The  correlation  Length,  which  will  be  used  later,  can  easily  be 
shown  to  equal 


In  2 


(2.3) 


for  the  Gaussian  curve  $  corresponding  to  $  .  This  formula 

n  n 

enables  us  to  obtain  a  good  estimate  for  the  correlation  length 

of  i  .  We  conclude: 
n 

The  correlation  length  of  the  cardinal  spline  base  function 

of  degree  n  is  approximately  proportional  to  vn+1  . 

Another  evident  fact  concerns  the  integral  of  <i>  :  we 

n 

recall  that  ;  =  y  ,  can  be  considered  a  density  function; 

n  n  +  1  J 

the  integral  is  the  distribution  function  and,  consequently, 


<t>  (x)dx  =  1  .  (2.9) 

j  n 

—  cc 

How  does  a  behave  if  n  goes  to  infinity?  We  conclude  from 
n 

(2.7)  that  the  maximum  of  j  (which  is,  of  course,  attained 

n  _ 

at  the  origin  x=0  )  decreases  like  l//n+l  .  Due  to  (2.8)  the 
product 

hi  (0)  =  /ini.  =  const.  (2.10) 

bet  us  denote  the  limit  function  by  j (x)  , 


5  (x) 


lim 


;  6 

■/  tt  (n+1 ) 


6  x2 
n+  1 


n-*-o> 


e 


(2.11) 
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Because  of  (2.8)  and  (2.10)  4 (x)  is  a  kind  of  layer  on  the 

real  line  with  infinitesimal  function  value  and  unlimited  length 
but  such  that  its  volume  equals  1  , 

00 

^  T (x)dx  =  1  (2 .12) 


Since  all  4>n's  are  frequency  functions,  the  convolution  between 

a  function  f  and  $  results  in  a  smoothed  function  f  , 

n  n 


f  (x)  =  f(x)  *  $  (x) 
n  n 

00 

=  1  f (x‘ ) $  (x1  -  x)dx’ 

J  n 


If  n  goes  to  infinity,  it  is  obvious  that 


'  f (x* )  5 (x'  -  x)  dx'  =  f 

J 

—  oo 


(2.13) 


where  7  denotes  the  mean  value  of  f  .  We  conclude: 

The  result  of  a  convolution  between  a  function  f  and 
the  distribution  <5  equals  the  mean  value  of  f  . 

This  distribution  will  be  used  later. 


SPECTRAL  PROPERTIES  OF  GAUSSIAN  SPLINE  APPROXIMATION 


The  Fourier  transform  of  i>  ,  the  Gaussian  curve  aoprox- 

n 

.mating  the  cardinal  spline  base  function  of  degree  n  (equu- 
:ion  (2.7))  ,  is  given  by  (Gradshteyn,  1971,  p.  494,  No.  3.896) 


>  (t) 

n 


/  6 
^  ( n+ 1 ) 


e  e 


dx  , 


5  (■>) 
n 


(2.14) 


Equation  (2.14)  can  also  be  written  as 


o  n  1 1 
(  T  1  )2 

i  { n )  =  !  e  6  ] 


2.14') 


and  shows  once  more  the  fact  that  t  resembles  the  n-fold 

n 

convolution  of  Gaussian  curves  , 


■t1  =  (  $n)  ,  $  o  (  n )  =  e 


(2.14") 


(Recall  the  corresponding  spectrum  of  the  base  spline  (subsection 
1.5.15), 


Jn  =  ' n+1  $0 (n) 


sinii 


(2.15) 


The  spectra  ♦  and  «  are  shown,  for  various  degrees  n  ,  in 
Fig.  2.2,  an  upper  bound  for  their  absolute  differences  in  Fig.  2.3 


i - 1 - 1 - — i - 1 - r 


i - 1 - r 


i - 1 


! . 50  1.00  1.50  2.00  2.50  3.00  3.30  4.00 


FIG.  2.2  Fourier  transforms  of  cardinal  spline  base  functions  i 
and  corresponding  Gaussian  functions  t 
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10  15  20 


degree  n 


25 


FIG.  2.3  Upper  bound  max  <t>  (n)-J>  (n) 


From  equations  (2.14)  and  (2.15)  as  well  as  from  Fig.  2.2  and 

2.3  we  conclude  that  both  spectra  4  and  4  become  identical 

n  n 

if  n  goes  to  infinity 

lim  i  =  lim  4 
n  n 

n-*'1’  n-*°“ 

The  limit  spectrum  equals  1  for  the  zero  frequency  and  vanishes 
for  all  other  frequencies.  Note  that  has  no  volume  (mass), 

therefore  it  is  not  an  impulse.  We  call  this  "function"  D ( n )  : 

D ( n )  :  =  lim  4  ( n )  .  (2.16) 

n 

n+°° 

Recalling  the  definition  of  the  distribution  6 (x)  (equation 
(2.11)),  it  is  obvious  that  D  is  the  Fourier  transform  of  5  , 

oo 

D(n)  =  [  T(x)  e  "i2"rixdx  .  (2.17) 

—  ao 

The  inverse  transform  of  D  is  the  distribution  5  , 

TO 

6  (x)  =  I  D(n)  e  i2TTnxdn  .  (2.18) 

J 

—  oo 

Therefore,  T(x)  and  D(n)  can  be  considered  as  a  Fourier  trans¬ 
form  pair. 

The  distribution  T  is  very  closely  related  to  the  Dirac 
distribution  -  as  summarized  in  the  following  table. 


I 
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The  mutual  relation  between  both  distributions  becomes  particularly 
clear  by  the  relation 


30  00 

r  _  r  __ 

1  D ( o ) 6 ( n ) d n  -  '  D(n)o(n)dn  =  1  ,  (2.19) 

;  j 


which  can  be  verified  using  the  above  relations.  Both  distributions 
can  be  considered  as  extreme  cases  of  kernels  of  a  moving  average 
integral  operator  M  , 


JO 

r 

M  (■)=:(•)  '5  dx  (  2.20a) 

n  J  n 

—  CO 


i 


with 


e.g  .  ,  defined  by 


(x)  =  —  e  ;  (2.20b) 

n  v  *t  n 

if  n  goes  to  zero,  we  obtain  the  Dirac  distribution  d  ( x. )  ,  if 
n  goes  to  infinity,  we  obtain  the  5 (x)  -  distribution;  both 
limit  operators  Mc  and  M  degenerate  into  functionals;  the 
evaluation  functional  M0  (which  "averages"  the  function  with¬ 
in  an  infinitesimally  small  interval) ,  and  the  averaging  functional 
Mm  (which  averages  the  function  on  an  infinitely  large  interval) . 
It  is  due  to  its  relatively  unimportant  properties  that  the 
distribution  5 (x)  is  not  considered  in  the  literature. 

Let  us  now  consider  the  Fourier  transform  of  the  Gaussiar. 
cardinal  (sampling)  function  corresponding  to  the  cardinal  (samp¬ 
ling)  function  of  degree  n  , 


v„(x)  =  i  *n<x  “  k>  • 


(2.21) 


According  to  subsection  1.5.1,  the  Fourier  transform  Y  (n)  is 

n 

given  by 


’  ( ■' )  =  P  ( n )  }  a  cos2-rkn  , 
c  n  k 


(2.22) 


whicn,  in  turn,  can  simply  be  transformed  to 


$  n  (  n  ) 
Y  (n)  =  - 

n  a  (n) 

n 


(2.23a) 


UPPtH  BOUND 
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with  $  (n)  given  by  equation  (2.14)  and  a_  ( r, )  by 


a  ( n ) 
n 


l  *  (k)cos2vkn 


(2.23b) 


Explicitly  written  f  (n)  is  as  follows 


~  (  tt  n ) 


(  r.  +  1  ) 


■f  ( n ) 
n 


JO  - 


6  k  ‘ 


/  iTnU)  +  2ie  n+‘  cos2irkn 

!.  1 


(2.24) 


The  following  graph  shows  the  maximum  absolute  differences  between 

r  and  f  for  different  degrees  n  . 
n  ,  n 
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FIG.  2.4  Upper  bound  max  |l'n(n)-j'n(n) 

Vn 


How  does  the  Limit  spectrum  •?  (a)  look  like?  Figure  2.5 

‘•no  *' 

demonstrates  the  behavior  of  ?  for  large  n  (n  =  10,  20,  ... 

n 

obviously  v  is  a  rectangular  pulse 

1  for  j  n  ;  < 

■v  (r)  =.  <  1 2.25 ) 

TO  1 

0  else 

Comparing  v  with  (equation  (1.64)),  we  see  that  both 

are  the  Fourier  transform  of  the  function  sinirx/iTX  , 


T  (  n  )  =  f  (  a  ) 


co 

e 

i 


sin^x 

7TX 


-  i  2  tt  n  x 
e 


dx 


(2.26a) 


and ,  consequently , 


y  (x)  =  U  (x) 


sinTTX 


~  x 


(2.26b) 


This  result  has  already  been  obtained  by  a  different  approach  in 
section  1.5. 


'f  ( n  )  ,  n 


FIG.  2.5  Convergence  of 
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2 . 2  DEVIATION  OF  GAUSSIAN  CARDINAL  BASE  AND  SAMPLING  CURVES 
FROM  THE  CORRESPONDING  SPLINES 


It  is  of  considerable  interest  to  know  a)  how  much  the 
cardinal  spline  base  functions  differ  from  the  corresponding 
Gaussian  curves,  and  b)  how  much  the  cardinal  sampling  splines 
differ  from  the  corresponding  Gaussian  ones.  The  first  question 
has  a  simple  answer.  Generalizing  equation  (1.68)  for  real  argu¬ 
ments,  one  arrives  at  the  following  expression  for  the  cardinal 
spline  base  function  of  degree  n  , 


9  ( X) 


.  ,'n+l  i  ■  N 

int!—  -t-  I  x  i  j 


n : 


(-1) 


j=0 


i'n+1)  'n+1  .In 

1  j  ;•  2 


for  0< 


x 


;  0  else  ; 

l 


5  (x)  is  given  by  (2.7)  , 


(2.27) 


l  ( x ) 


it  ( n+ 1 ) 


e 


fc  x^~ 

n+  1 


which  makes  the  calculation  of  i  -  p  an  easy  task.  The  graph 

n  n 

in  Fig.  2.6  illustrates  the  convergence  behavior  of  the  deviation 


5  (x)  -  ■>  (x) 
n  n 

max  -  ,  n  -<■  «  . 

/x  •;  (O) 

n 


Another  estimate  of  an  upper  bound,  which  is,  of  course,  more 
pessimistic,  can  be  obtained  via  the  spectra, 
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H  + 
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DEGREE 


The  second  question  of  how  much  the  cardinal  sampling  splines 
differ  from  the  corresponding  Gaussian  curves,  can  virtually  only 
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be  answered  by  employing  the  ^pectra.  The  absolute  difference, 
as  above,  can  be  represented  in  the  following  way 


v  (x)  -  '1  (x) 
n  n 


( r, )  -  r  _  (  n )  cos  2  -  ",xdn 


(2.29, 


2  \'t  (n) 

i  '  n 


V  (n)  dr 
n 


with  i  (-)  =  ■>  (n)/a(n)  and  r  ( •> )  =  1  (0/a(n)  ;  the  exolic- 
n  n  nr. 

it  expressions  are  given  by  (1.71)'  and  (2.24).  The  graph  in 
Fig.  2.7  demonstrates  the  convergence  behavior  of 

max  : vn (x)  -  vn (x) ;  ,  n  ; 

v  X 

(note  that  ;  (0)  =  1  for  all  n  ;  :  (0)  ,  however,  is  less 

n  n 

than  1  for  n  *  1  .) 

Fig.  2.8  shows  cardinal  sampling  splines  of  various  degree 
and  the  corresponding  Gaussian  sampling  curves. 


.  loo  -I 


aerjree 


FIG.  2.7  Convergence  of  max  _  <  x )  —  , 
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For  prediction  problems  in  geodesy  the  Gaussian  covariance 
function  is  considered  not  as  important  as  the  covariance  function 
of  Hirvonen  type  which  has  the  following  general  structure 

-  v  2-~a 

$(:<)  =  L1  +  (ij)  _  (2.20) 

Three  values  for  a  are  of  particular  interest  : a  =  1 / 2 ,  1,  3/2; 
for  a  detailed  discussion  see  (Moritz,  1976)  .  The  parameter  a  is 
related  to  the  correlation  length  h  through 


Fcr  a  cardinal  data  distribution  the  corresponding  sampling 
function  y(x)  is  obtained  in  the  same  way  as  discribed  for 
the  Gaussian  case.  For  the  sake  of  completeness  the  sampling 
functions  of  Hirvonen  type,  dependent  upon  the  correlation  length 
h  ,  are  presented  for  a  =  1/2,  1,  and  3/2;  the  correlation 
length  has  been  selected  as  h  *  0.1,  1,  and  10  in  order  to 
point  out  the  essential  features.  Comparing  Fig.  2.9  (Hirvonen) 
with  Fig.  2.8  (Gaufi)  and  taking  into  account  equation  2.8,  we 
conclude  ;  to  each  Hirvonen  sampling  function  an  almost  indis¬ 
tinguishable  Gaussian  sampling  function  exists;  the  correlation 
length  of  the  Gaussian  base  function,  however,  has  to  be  signifi¬ 
cantly  smaller  than  the  one  of  the  corresponding  Hirvonen  base 
function.  This  statement  is  valid  at  least  for  the  Hirvonen 
functions  considered  here  (ot  =  1/2,  1,  3/2).  This  fact  has  a 
very  strong  impact  on  interpolation  error  estimation  discussed 
m  chapter  5  .  With  increasing  correlation  length  h  the 
Hirvonen  sampling  functions  approach  the  sinirx/ix  function; 
it  shares  this  feature  with  the  Gaussian  sampling  function. 
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3 .  SPLINE  INTERPOLATION  VERSUS  LEAST-SQUARES  INTERPOLATION 

In  the  foregoing  chapters  it  has  been  shown  how  splines 
of  different  degree  are  related  to  corresponding  Gaussian  curves. 
The  link  is  provided  by  the  central  limit  theorem  or  by  the 
theory  of  linear  systems  in  cascade  (Papoulis,  1968,  p.  50) . 

All  spline  interpolations  on  a  homogeneous  set  of  regularly  dis¬ 
tributed  data  along  the  real  line  lead  to  the  following  form 


f  (x)  =  l  f  v  (x  -  k)  ;  (3.1) 

n  u  k  n 

—  CO 

(for  the  sake  of  simplicity  a  data  interval  equal  to  unity  has 
been  assumed.)  It  is  essential  to  point  out  that  the  function 
is  a  sampling  function  whose  support  is,  in  general,  the 
whole  real  line  (if  one  ignores  the  countably  infinite  number 
of  zeroes);  only  l>  0  and  ^ .  have  a  limited  support  of  1 
and  2  ,  respectively.  -j;  is  ( n- 1 )  times  continuously  differ¬ 
entiable  and  consists  of  piecewise  polynomials  of  degree  n  . 

is,  apart  from  the  cases  n=0  and  n=l  ,  not  directly 
given;  it  is  rather  the  result  of  a  linear  transformation  of 
originally  defined  base  functions  t>_(x)  such  that 


r  ( x )  = 


t .  t  ( x  -  k  ) 
k  r4 


The  coefficients  ■  <.  can  be  determined  from  the  sampling  proper 

K 

ties  of  '1/  .  This  leads  to  the  solution  of  tne  linear  system 

n 


The  system  can  be  solved  by  means  of  frequency  domain  methods 
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and  provides  simple  solutions  for  •  t  •  if  rr-3  .  Since  the  in¬ 
finite  matrix  v  :  =  •' t  (7  -  k) }  is  of  Toeplitz  type,  its  inver¬ 
se  i>~;  is  also  a  Toeplitz  matrix;  consequently,  '.he  infinite 
solution  vector  1  =  /jlk  represents  any  row  or  column  of 
: " *  is  completely  defined  by  a  , 

:t  =  t  ~  -  e 

The  interpolation  function  f  (x)  can  also  be  formulated  as  a 

n 

linear  combination  of  the  base  functions  $  (x)  , 

n 

a  00 

f  (x)  =  l  6.  0  (x  -  k)  .  (3.2) 

n  <  n 

—  s> 

According  to  the  interpolation  property  there  must  hold 

SMJ)  -  f  -  I  -  K) 

■“  00 

which  yields  the  solution  vector  6 

fi  =  *> "  1  f  ; 

is  given  by  a  ,  therefore  3  is  known  (f  is  the  data 
vector) ,  and  we  obtain  as  solution  to  the  interpolation  problem 

f  (x)  =  t  (x-k){i  ( k  -  1 )  1  -  f , 

r.  n  n  x. 

(The  summation  over  k  and  1  is  self-evident.)  Briefly, 


In  order  to  be  absolutely  clear,  let  us  once  more  explain,  what 


o  o 


these  symbols  mean:  f  is  the  data  vector,  f 

function  value  at  the  point  P  (x=P)  ,  ?p  is 

lated  base  functions  {9  (P  -  k) }  evaluated  at 

n 

;  is  the  matrix  of  transferred  base  functions 
at  the  points  x=  1  ,  1,  k  =  -1,  0,  1,  . 


the  interpolated 

the  vector  of  trar.s- 

the  point  P  x=P)  , 

9  ( 1  -  k )  evaluated 

■  n 


Formally,  the  spline  function  interpolation  solution  (3.3) 
is  identical  with  the  least-squares  interpolation  solution  (e.g. 
Moritz,  1980,  p.  81) 


f  =  cV:f  .  (3.4) 

If  we  consider  the  covariance  function  C(x)  a  "base  function", 
the  relation  becomes  evident;  all  what  has  been  said  above  about 
the  base  functions  >|>  (x)  can  be  literally  copied  for  the  covari¬ 
ance  (base)  function  C(x)  : 

The  interpolation  function  fc(x)  can  be  expressed  as  a 
linear  combination  of  sampling  functions 


f(x)=/fv(x-k)  ,  (3.5) 

C  w  K  C 


or  by  a  direct  linear  combination  of  the  covariance  (base)  functions 


f  (x)  =  7  3,  C(x  -  k)  (3.6) 

c  k 

— •no 


The  interpolation  property  leads  to’  the  solution  of  the  linear 
system 


f 

c 


( j) 


f  = 


C(j 


k) 


which  yields  the  solution 


=  =  G'-f  ; 

the  interpolation  function.  f  (x)  is  finally  given  by 

f(x>=C(x-k)'Cfk-l)}-if 

Note  that  the  least-squares  sampling  function  ;^(x)  is  expressed 
as  a  linear  combination  of  the  covariance  functions;  the  coeffi¬ 
cients  of  the  linear  combination  are  the  elements  of  the  inverse 
covariance  matrix  C~  !  , 

„•  (x  -  k)  =  C*!c(x  -  j)  . 

c  kj 

This  can  immediately  be  verified  by  equating  (3.5)  and  (3.6),  ob¬ 
serving  =  C~ 1 f 

The  two  solutions  share  many  common  features  (also  the  co- 
variance  matrix  is  of  Toeplitz  type,  etc.);  there  is  one  very 
essential  difference,  however,  between  both  solution,  which  is 
due  to  the  structure  of  the  base  function  :  the  covariance  (base) 
function  has  an  unlimited  support,  therefore,  the  covariance  ma¬ 
trix  is  a  full  matrix;  the  spline  base  functions  have  all  a  limitei 
support,  the  corresponding  spline  matrices  are,  therefore,  band 
matrices  with  a  bandwidth  N  depending  on  the  degree  n  of  the 
spline , 


N  =  2  •  int  (~>  +  1  (3.7) 

l. 

Band-structured  matrices  can  be  inverted  much  faster  than  full 
ones  using  a  modified  Cholesky  decomposition  algorithm;  needless 
to  say,  the  speed  depends  on  the  actual  dimension  of  the  matrix 
and  on  its  bandwidth  :  if  M  denotes  the  dimension  of  the  matrix 


ar.d  n  the  degree  of  the  spline,  the  inversion  time  is  proportion¬ 
al  to  int(^)  ,  •  '  M  as  compared  to  M‘  for  a  full  matrix 

Forsythe  and  Holer,  1967);  with  other  words,  a  band-matrix  inver¬ 
sion  is  on  the  order  of  M/int(^)  faster  than  the  inversion, 
c:  a  full  matrix.  This  fact,  among  others,  makes  spline  interpo¬ 
lation  particularly  attractive. 

Let  us  put  a  very  crucial  question;  Can  a  spline  interpo¬ 
lation  replace  a  least-squares  interpolation? 

In  chapter  2  it  was  demonstrated  how  splines  of  increasing 
degree  approach  the  Gaussian  curve.  If,  therefore,  the  covariance 
function  is  of  Gaussian  type,  and  if  a  spline  of  sufficiently  high 
degree  is  chosen,  we  can  expect  that  least-squares  interpolation 
and  spline  interpolation  give  practically  the  same  result;  this 
should  be  true  not  only  for  the  interpolated  value,  but  also  rcr 
the  predicted  interpolation  error.  The  degree  n  of  the  spline 
depends  only  on  the  correlation  length  of  the  corresponding  Gaus¬ 
sian  covariance  function  through  equation  (2.8)  , 


n 


int 


1 


(3.8) 


(In  this  context,  h  refers  to  a  data  spacing  equal  to  unity.) 
Since  the  correlation  distance  of  spline  base  functions  increases 
only  with  the  square  root  of  the  degree  n  ,  in  general,  a  very 
high  degree  spline  has  to  be  chosen;  as  a  consequence ,  the  band¬ 
width  of  the  spline-approximated  covariance  matrix  will  increase: 
according  to  equation  (3.7)  the  bandwidth  increases  linearly  with 
the  degree  n  and,  considering  (3.8)  ,  quadratically  with  the 

normalized  correlation  length  of  the  covariance  function. 

E.g.  :  Assume  1000  data  (M  =  1000)  ;  a  full  inversion 

on  a  very  fast  AMDAHL  470  system  would  require  some  20  minutes 
of  CPU-time  ;  the  following  estimates  are  obtained  for  correspond¬ 
ing  spline  solutions: 


h 


n 


I  CPU  (sec) 


1 

7 

0.01 

2 

33 

0.3 

i 

3 

77 

1  .6 

4 

!  137 

i  5.3 

5 

;  215 

j  12.8 

The  band  matrix  becomes  a  full  one  if 


=  10.8  in  our  example.  If  the  degree  n 
increases,  the  calculation  of  a  spline  function  value  according 
to  (1.68)  becomes  expensive.  This  should  not  cause  any  problems, 
however,  since  in  practical  applications  one  will  generate  a 
high  density  table  for  the  base  function  once  for  all  and  find 
an  actual  function  value  by  a  simple  table  interpolation  function; 
this  is  also  a  common  practice  in  applications  of  least-squares 
collocation,  where  covariances  are  usually  tabulated;  an  inexpen¬ 
sive  interpolation  algorithm  is  described  in  (Siinkel,  1973,  1979)  . 

So  far  we  have  assumed  uniform  data  distribution  which  is 
hardly  ever  available  in  practical  applications.  Spline  as  much 
as  least-squares  interpolation,  however,  is  by  no  means  restricted 
to  regular  distributions.  For  irregularly  distributed  data  the 
basic  equations  (3.3)  and  (3.4),  respectively,  still  hold,  although 
the  basic  features  of  the  spline-  and  corresponding  covariance 
matrices  disappear  :  the  matrices  are  no  longer  of  Toeplitz  type 
(C  ^  C(j  -  k)  ,  in  general);  therefore,  the  spline  matrix  is 
no  ionaer  a  band-matrix;  it  is,  however,  still  a  sparse  matrix 
with  the  number  of  zeroes  approximately  equal  to  that  in  the  cor- 


M  +  l  ,  , 
h  .  - - -  ln2 

/  6 


which  corresponds  to  h 
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responding  band-matrix  (the  correspondence  referring  to  an  average 
data  spacing ) . 

If  a  least-squares  interpolation  is  replaced  by  a  corre¬ 
sponding  spline  interpolation,  it  is  required  that  the  spline 
base  function  (which  replaces  the  covariance  function)  be  positive 
semi-definite  which,  in  terms  of  the  spectrum,  means  that  the 
spectrum  of  the  introduced  base  function  has  to  be  non- negative . 

We  know  from  chapter  1  that  the  spectrum  of  *  is  given  by 

¥  n 


i 


i 


r  '  n*  1 

i  sm-  n 

I  7  rj  i 


since  simn/n  is  positive  and  negative,  we  need  an  even  power 
in  order  to  keep  the  spectrum  non-negative;  an  even  power  of 
sinirn/-rn  corresponds  to  an  odd  degree  spline  base  function.  There¬ 
fore,  we  conclude  that  only  odd  degree  spline  iterpolations  can 
be  used  to  replace  a  corresponding  least-squares  interpolation. 

An  interesting  question  arises  in  connection  with  the 
above  described  relation  between  spline  and  least-squares  inter¬ 
polation  :  It  can  be  shown  that  the  spline  interpolant  of  degree 
2n-l  minimizes  the  quadratic  functional 

co 

’ f ! " | i 2dx  =  minimum  (3.9) 

j  .  X  ) 

—  JO 

with  f' "  denoting  the 

( x ) 

et  al .  ,  1967)  .  Isn't  it 

problem  of  the  spline  to 
squares  interpolation, 

f TC“ ! f  =  minimum  ?  (3.10) 


n'th  order  derivative  of  f(x)  (Ahlberg 
possible  to  relate  this  norm  minimization 
the  norm  minimization  problem  of  least- 


i 


I 


Let  us  express  f  in  equation  (3.9)  by  the  spline  given  in 
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equation  (3.3)  , 

f  ( x )  =  :> 7  ( x )  ->  ‘  *  f  , 

f(x)  =  :,,7(x)  <t>- ‘  f  ;  (3.1i: 

with  (3.11)  the  minimum  condition  (3.9)  is  transformed  into 


x 


—  oo 


f 

3 


-1 


OO 

r 


(ni 


k)* (n}  (X 


1 )  dx 


A  * 

"lir 


i. 

ra 


(3.12) 


Here  we  have  explicitly  expressed  the  dependence  of  j  on  the 
difference  x-k  and  x-1  ,  respectively.  With  (3.12)  the  minimum 
condition  (3.9)  can  be  transformed  such  that  its  dependence  cn 
the  data  is  obvious, 


fT-t>~  1  Ke‘  1  f  =  minimum  , 

00 

K,  :  =  i{ri\x  -  k )  <$  1  n '  ( x  -  1 )  dx 
x  1 


(3.12)  • 


It  is  clear  that  (3.12)  ’  would  reduce  to  the  least-squares  inter¬ 
polation  condition  of  the  form  (3.10)  if  the  matrix  K  would 
equal  3  .  Does  it? 

K  represents  obviously  a  convolution;  therefore,  the 
elements  of  K  depend  only  on  the  difference  k-1  , 

K.  =  K(k  -  1)  ; 

x  1 

consequently,  K  is  a  Toeplitz  matrix.  The  function  3{n)(x) 
is  the  n-fold  derivative  of  the  spline  base  function  $ „  (x)  - 
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in  the  notations  above  we  have  consistently  suppressed  the  degree 
2n  -  1  -  which  has  a  support  of  length  equal  to  2n  ;  since 
:  (  1  x '  +  e )  <  .?  (  x  )  ,  e  >0  1  1 )  within  the  support  x  r.  , 

it  follows  that  the  support  of  the  derivative  py"‘  is  also 
equal  to  2n  .  The  convolution  p ( n '  *  has  twice  the  support 

of  ?  ;  consequently  K  cannot  equal  i  and  the  answer  of  the 

above  question  is  in  the  negative. 

Furthermore,  utilizing  (1.68)  ,  it  can  be  shown  that  the 

sum  of  each  row  (or  column)  of  K  equals  zero,  whereas  the  sum 
of  each  row  (or  column)  of  $  equals  1  .  Such  kind  of  matrices 
are  known  as  "stochastic  matrices"  (Zurmiihl,  1964  ,  pp.  221-224)  . 

Such  kind  of  matrices  are  known  to  have  a  simple  eigenvalue 
which  equals  the  row-sum;  the  corresponding  inverses  have  an  eigen¬ 
value  which  equals  the  reciprocal  value  of  the  row-sum.  The  product 
of  two  stochastic  matrices  has  an  eigenvalue  which  equals  the 
product  of  their  row-sums.  Therefore,  >  and  have  an  eigen¬ 

value  equal  to  1  ,  K  has  a  vanishing  eigenvalue  and,  consequently, 

the  product  has  also  a  vanishing  eigenvalue.  With  other 

words,  the  "metric"  is  singular.  As  a  consequence,  there 

are  data  vectors  f  which  are  not  identically  zero  and  have  still 
a  vanishing  norm  in  the  space  defined  by  the  above  metric.  Such 
norms  are  called  "semi-norms";  they  are  characteristic  for  spline 
functions.  E.G.:  The  norm  minimization  problem  (3.9)  with  n=l 
(linear  spline)  is  "blind"  with  respect  to  a  constant  function; 
the  norm  minimization  problem  with  n=2  (cubic  spline)  is  "blind" 
with  respect  to  a  constant  function  and  with  respect  to  a  linear 
function  with  constant  derivative;  in  general,  the  norm  minimizaticn 
problem  (3.9)  for  the  spline  of  degree  2n-l  is  "blind"  with 
respect  to  all  functions  g(x)  and  all  linear  combinations  thereof, 
having  constant  derivatives  D  j ( x )  ,3=0,  ...,  n-1  ;  these  are 

obviously  polynomials  of  degree  n-1  .  This  "blindness",  however, 
does  not  have  any  consequences  in  the  case  studied  here,  as  will 
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be  argued  ir.  the  sequel: 

The  minin'.! nation  problem  (3.9)  is  admissible  in  a  space 
of  functions  raving  square  integrable  n’th  order  derivatives, 
called  I,‘.‘  ;  then  it  follows  that  the  data  vector  has  to  be  an 

element  of  11  ,  the  space  of  square  integrable  n 1 th  order 

divided  differences.  (For  an  explanation  of  divided  differences 
see,  e.g.  Davis,  1975,  pp.  40,  64  ff . . )  It  can  be  shown  that 
under  these  presumptions  and  the  interpolation  condition,  there 
exists  a  unique  solution  to  the  minimization  problem  (3.9)  such 
that  tile  interpo lating  function  f(x)  is  a  spline  of  degree 
2n-i  (Schoenberg,  1973  ,  p.  59 f f . ) 

The  norm  for  the  least-squares  solution  differs  consider¬ 
ably  from  the  norm  (3.9)  :  since  the  covariance  matrix  is  posi¬ 

tive  definite,  the  norm  (3.10)  vanishes  if  and  only  if  the  data 
vector  vanishes  identically.  In  chapter  2  it  was  shown  that  the 
spline  solution  approaches  the  least-squares  solution  if  the  cor¬ 
relation  length  (and,  accordingly,  the  degree  of  the  spline)  in¬ 
creases  .  How  can  two  thus  different  minimization  problems  lead 
to  practically  the  same  result? 

The  key  is  obviously  given  by  the  spline  norm  on  the  one 
hand  and  the  relation  between  splines  and  the  Gaussian  function 
on  the  other  hand  :  the  spline  of  degree  2r.-l  minimizes  the 
integral  of  its  squared  n 1 th  derivatives  and  interpolates  the 
data.  As  stated  before,  this  seminorm  minimization  is  blind 
with  respect  to  a  polynomial  of  degree  n-1  ;  due  to  the  inter¬ 
polation  condition,  however,  this  polynomial  is,  so  to  say, 
"smuggled  in";  the  higher  the  spline’s  degree,  the  less  important 
is  the  corresponding  norm  minimization  and  the  more  pronounced  is 
the  (single)  polynomial's  influence  on  the  behavior  of  the  inter¬ 
polation  function.  (It  should  be  evident  that  the  minimization 
of  the  squared,  e.g.  5 ' th  order,  derivatives  contributes  only 
very  little  to  the  shape  of  the  function.)  In  the  limit,  if  the 
spline's  degree  goes  to  infinity,  its  corresponding  samplina 
function  equals  the  interpolation  function 
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f  (x) 
» 


sin^i (x-k) 
t (x-k) 


(3.13) 


which  has  been  shown  in  section  1.5.  Therefore,  if  our  conclusions 
concerning  the  polynomial-like  behavior  of  high  degree  splines 
are  valid,  the  above  spline  interpolation  of  infinite  degree  snould 
behave  like  a  single  polynomial  of  infinite  degree  interpolating 
the  cardinal  data  vector.  More  specifically,  the  polynomial  inter¬ 
polation  needs  to  be  expressible  as  a  linear  combination 


co 


with  y  Jx)  being  a  sampling  polynomial, 


(3.14a) 


!  1 

at 

X 

(X)  = 

i  0 

1 

at 

X 

¥■  o 

else 

(3.14b) 


This  sampling  polynomial  is  a  Lagrange  polynomial  of  infinite 
degree  given  by  (e.g.  Moritz,  1978,  p.  3)  , 


.  .  ,  _  ...  (x-k+2)  (x-k+1 )  (x-k-1)  (x-k- 2)  ... 

-  ~  . . .  2  •  1  •  (-1)  •  (-2;  ... 

(3.15) 

Whittaker  (1935,  p.  62ff.)  verified  that,  in  fact,  this  Lagrange 
sampling  polynomial  equals  the  simi/v  -  sampling  function, 


3inu< 
-  x 


(X) 


(3.16) 
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This  rer.arkablo  relation  provides  'is  with  tne  perfect  proof  that 
high  degree  splines  tend  to  behave  like  single  polynomials;  the 
spline  of  infinite  degree  is,  in  fact,  n  single  polynomial  of  . 
infinite  degree  —  the  minimisation  of 

r 

1  r  c  1 .  n  )  .  -  .  . 

lira  f,  * dx  =  minimum 

i  (  X  ! 

n 

—  or? 

looses  its  meaning  completely. 

Let  us  return  to  the  above  question.  It  was  shown  in 
chapter  2  that  splines  approach  the  Gauss : an  curve  with  increasin 
correlation  length  and  correspondingly  increasing  degree.  Since 
both  interpolations,  the  spline  and  the  least-squares  interpolation 
have  the  same  structure  (equations  (3.3)  and  (3.4),  we  conclude 
that  the  norm  minimization  (3.10)  of  the  least-squares  solution 
becomes  meanlnaless  if  the  correlation  length  increases  in  very 
much  the  same  way  as  the  norm  minimization  (3.9)  of  the  spline 
solution  becomes  meaningless  with  increasing  degree.  Therefore, 
both  solutions  have  to  behave  in  the  same  way  in  this  limit  case  - 
like  a  polynomial  interpolation. 
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4 .  8JERHAMMAR  INTERPOLATION 


Both  the  least-squares  interpolation  and  spline  interpola¬ 
tions  share  a  common  feature:  the  determination  of  the  interpola¬ 
tion  function  requires  an  inversion  of  a  matrix;  its  size  equals 
the  number  of  data  for  the  case  of  least-squares  interpolation, 
it  reduces  to  a  band  matrix  for  the  case  of  spline  interpolations 
with  a  bandwidth  depending  on  the  degree  of  the  spline  (or  on  the 
correlation  length).  With  other  words,  the  sampling  functions  are 
obtained  through  an  inversion  process ,  a  usually  expensive  under¬ 
taking.  Is  there  a  method  which  avoids  this  inversion  —  can't 
the  sampling  functions  be  replaced  by  an  a  priori  weight  function? 
This  question  has  been  considered  by  various  authors  like  Bjerham- 
mar  (1973)  or  Shepard  (1964) .  We  will  consider  in  detail  only  the 
method  proposed  by  Bjerhammar  because  it  was  and  is  still  being 
considered  a  particularly  useful  tool  for  data  interpolation  by 
members  of  the  geodetic  community. 

Bjerhammar  proposes  to  define  the  interpolation  function 
as  a  weighted  average  of  the  data  with  weights  depending  on  the 
distances  between  the  prediction  point  and  the  data  point, 


-  M 

f  ,  =  '  f  :i  (  '  X-X  '  )  . 

••an  ra 

m  -  1 


(4.1) 


in  this  context  f  denotes  the  weight  function;  it  replaces  tne 
a amp  line  function  used  before.  In  the  case  studied  here,  an  in¬ 
finite  data  set  is  given,  defined  at  the  cardinal  numbers, 


f  ( x)  =  y  f  (  x-m ;  )  . 

m 


(4.1)' 


7  9 


Since  fix) 
that  (x) 
i unction , 


interpolates  the  data  sequence  f  >  ,  it  follows 

m 

needs  to  fulfil  the  conditions  of  a  sampling  (weight) 


?(j>  = 


(4.2) 


8 ( ! x-m ! )  =  1 


The  weight  functions 


6  (x)  : 

P 

are  proposed  by  Bjerhammar*;  the  positive  quantity  p  will  be  de¬ 
noted  power  of  prediction. 

Particularly  simple  is  the  Bjerhammar  sampling  function 
if  p  =  1  : 


l  (x-k)“(pM) 


3;  (X) 


1 


no 


—  CO 


(x-k)“2 


(4.3) 


The  infinite  sum  in  the  denominator  can  be  expressed  by  a  closed 
expression  (Gradshteyn,  1971,  p.  50,  Mo.  1.422), 


-  1  t  V- 

-•»  ( x-k )  •  isinnxi 


)  ...  Here  the  smoothing  quantity  has  been  neglected. 


80 


therefore,  Si(x)  is  given  by 


(x)  = 


i  sinirx  i 

\  TTX  J 


(4.3/' 


an  expression  which  is  already  familiar  to  us  :  £;(x)  is  nothing 

but  the  Fourier  transform  of  the  linear  interpolation  sampling 
function,  a  very  astonishing  result  (cf.  equation  (1.17))  , 


A  ■  ( x)  =  V  ,  (x) 


(4.4) 


(It  has  been  shown  in  section  1.2  that  V l  =  $ i  .) 

The  corresponding  interpolation  function  is  simply  given 
by 


Z ,  ,  r  ,  rsimr  (x-m)"  2 
f(x)  =  l  fm  ■  '■  ,(x-m)  J 

—  30  —  -* 


Since  the  sampling  function  s5  equals  the  Fourier  transform 
of  the  spline  sampling  function  v;  ,  it  follows  that  the  spectrum 
3  ■  of  r  i  equals  the  spline  sampling  function  'p  i  (confer 
section  1.2), 


B  (•,)  =  1  h  ■  (x)e_l27,r,xdx 


f 1  -  n ,  for  Ini  •  1 


else 


(4.5) 


B  .  (  n  )  -  v  (  n  ) 


(4.5)  ' 


conclude  : 


8  : 


v’e 

7’ir  5,'c 2"!«!ar  .*•:"?«  2  :,t  j  func  t  og’  vo  war 

"  ru  ’"o  •'?  r”  :  •'  cr  '  ir.a  .*  err  '  •'  >•  j  Fur:  >  •*  -:^r  >•->>: 

l* .  ;  v»  c?  , 

Let  us  now  consider  the  3jerhammar  sampling  function  of 
pewe  r  3  : 


1 

:’i(x)  - - 

-  00 

x"4  7  (x-k) 

—  oo 


(4.6) 


The  infinite  sum  in  the  denominator  can  be  expressed  by 


-  j  _  1  +  2  cos2  TX  [  TT  i  '* 

(x-k)"  ~  3  ^simrxi 

(Hansen,  1975,  p.  110,  No.  6.1.134);  using  the  relation 
2  cos  -  ttx  =  1  +  cos  2  rx  ,  6  3  reduces  to 


33  (x) 


_ 3 _  1  sin^rx)  '■* 

2  +  cos  2  n  x  v  -  x  j 


(4.6)  ' 


Tnis  is  exactly  the  Fourier  transform  'J  >  of  the  cubic  sampling 
spline  p  , 


3  ,  f  x) 


?  3  ( x )  . 


(4.7) 


I 
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The  corresponding  interpolation  function  is  given  by 


f  { x )  = 


sin- ( x-m) 


2  +  cos  2  ir(x-m)  i  it  (x-m 


'  ) 


Arguing  as  before,  the  Fourier  transform  B3  of  the  sampling  func¬ 
tion  33  equals  the  spline  sampling  function  y 3  , 

B3  ('i)  =  i|/j  (n)  .  (4.8) 


We  conclude: 

The  5,j  erhammar  sampling  function  of  power  .J  equals  the  Fourier 
transform  of  the  spline  sampling  function  of  degree  &  and  vice 
•-  rsa . 

How  does  the  Bjerhammar  interpolation  behave  if  the  power 
of  prediction  p  increases  beyond  any  limit? 

The  limit  sampling  function 


lim  s  (x) 

O 

p-^QO 


lim 

p->o> 


x 


p+  1 


or> 

l 


1 


x-k 


'1  p  +  l 

I 


can  be  simplified  to 


lim  S  (x) 

p>®  ^ 


I 

lim  — - - 

P-°"  7  1 

‘  • ,  kip+i 

-”  1  — 


I 
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r 

=  lim 


It  is  obvious  that  ^(x)  equals  1  for  |xl  <  0.5  ,  assumes 
the  value  0.5  for  x;  =0.5  ,  and  vanishes  for  Jx1  >0.5  , 

I  1  for  [  x  |  <  Vi 

3  (x)  =  i  '-'i  for  '  x  i  =  V'2  .  (4.9) 

I  0  else 
l 

This  is  exactly  the  step  function  introduced  at  the  beginning  of 
chapter  1  ,  it  is  at  the  same  time  the  Fourier  transform  of  the 
sampling  spline  of  infinite  degree.  Therefore,  we  conclude: 

The  behavior  of  the  Bjerhammar  interpolation  function  approaches 
zhat  of  a  step  function  "interpolation"  if  the  power  of  predictio' 
r  increases . 

Fig.  1.12  shows  Bjerhammar  sampling  functions  (odd  degree)  of 
various  degrees.  Fig.  4.1  is  an  example  of  Bjerhammar  interpola¬ 
tion  with  power  of  prediction  p  =  3  ;  the  data  are  located  at 
the  cardinal  numbers;  the  step-like  behavior  of  the  interpolation 
function  is  obvious. 


FIG.  4.1  Example  of  B jerhaminar  interpolation  (p  =  3) 

The  other  extreme  is  obtained  if  p  goes  to  -1  : 
tnis  "interpolation"  function  reproduces  (as  all  the  others)  the 
data;  in  between  the  data  it  assumes  a  constant  value  equal  to 
the  average  function  value  of  the  data;  it  is,  therefore,  a  con¬ 
stant  "function”  having  peaks  at  the  data  points. 

A  remarkable  and  important  property  of  Bjerhammar  interpo 
lation  is  the  following:  the  minima  and  maxima  of  the  interpola 
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tion  function  coincide  with  the  minimal  and  maximal  function 
values  of  the  data;  this  follows  immediately  from 


Summarizing  it  can  be  said  that  the  Bjerhammar  interpola¬ 
tion  with  odd  power  of  prediction  p  uses  the  spectrum  of  the 
corresponding  spline  sampling  function  as  its  sampling  function; 
these  functions  show  a  strong  tendency  to  approach  a  step  func¬ 
tion  with  increasing  p  ;  even  with  low  values  of  p  (cf.  Fig.  4 
a  step-like  interpolation  function  is  produced;  therefore,  it 
can  (if  at  all)  only  be  used  if  the  function  is  very  smooth,  if 
it  has  only  small  slopes,  and  if  it  is  well  sampled.  Bjerhammar 
and  others  recommend  a  power  p  =  2.5  ;  I  consider  this  value  as 
too  high  and  suggest  to  use  0.5  <  p  <  1.5  . 
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5.  INTERPOLATION  ERROR  ESTIMATES 


In  chapter  3  it  was  shown,  how  least-squares  sampling 
functions  change  with  the  correlation  length;  the  tendency  to 
approach  the  polynomial  behavior  corresponding  to  the  sin^x/'x 
sampling  function  was  obvious.  Particularly  interesting  is  the 
following  question:  how  good  does  the  data  set  represent  the 
function.  This  is  a  pure  sampling  problem  and  the  answer  can  be 
given  in  terms  of  the  least-squares  interpolation  error  estimate, 


mp 


=  C 


pp 


-  CTC-1C 
p  p 


(5.1) 


It  has  been  shown  that  the  least-squares  sampling  function  */(x)  , 

corresponding  to  the  covariance  (=  base)  function  <t>(x)  ,  equals 


'i>  (x-i)  =  i  <p  (i-  j  )  1  *  (x-  j  )  .  (5.2) 

With  this  notation,  the  predicted  error  variance  (5.1) 
reduces  for  an  infinite  data  set  at  the  cardinal  numbers  to 


nv  (x)  =  ;>  (0)  -  ?  (x-  j )  $  (x-  j )  .  (5.1)' 

j  =-°° 

The  interpolation  error  has  to  vanish  if  the  argument  x  is  any 
integer;  therefore,  in  this  case  the  infinite  sum  should  be  equal 
to  >  (O)  .  Observing  the  sampling  property  of  ,  ^  ( j )  =  3  , 

the  infinite  sum  reduces  to  a  single  component  ■*>  (O)  and 


m-  (k)  =  t (0) 
=  1  (0) 


l  ,  (k-j)  .■>  (k-j) 

j=-oo 

*(0)  =0 
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follows.  Due  to  symmetry,  the  maximum  interpolation  error  occurs 
just  in  between  the  data  (x  =  ■ '  k)  , 

M  :  =  max  m(x)  ’  *  mt1^  *  k)  .  (5.3) 

From  practical  considerations,  we  know  that  M  depends  on  the  cor¬ 
relation  length  h  of  the  underlying  covariance  function  :  a  small 
h  leads  to  large  interpolation  errors,  a  large  h  to  small  inter¬ 
polation  errors.  In  the  following, three  frequently  used  covariance 
models  are  considered  which  have  the  common  form 


»  (x) 


U 0)  I  1  + 


(5.4) 


Particularly  well-known  is  the  model  a  =  1  which  has  exclusively 
been  used  by  Hirvonen  as  early  as  1962.  The  other  two  models, 
x  =  'A  and  a  =  \  ,  have  been  discussed  by  Moritz  (1976)  because 
of  their  simple  relation  to  their  spatial  spherical  analogues. 

The  correlation  length  h  is  related  to  the  parameter  a  through 

./n*  1.4 

h  =  a  (2  -  I )  2 

which  reduces  to  h  =  a  for  the  Hirvonen.  model. 

The  following  Fig.  5.1  shows  the  maximum  interpolation  er¬ 
ror  M  'eq.  5.3) ,  depending  on  the  correlation  length  h  ,  for 
tnese  three  covariance  models;  a  variance  <t>(0)  -  1  has  been  as¬ 
sumed.  For  comparison  purposes  M  has  also  been  plotted  for  the 
corresponding  Gaussian  covariance  function. 


8b 


correlation  length  h 


FIG.  5.1  Maximum  interpolation  error  M 

The  following  conclusion  can  be  drawn  from  the  graphs: 
a)  the  maximum  interpolation  error  drops  fast  in  the  range  h <  2 
and  slowly  in  the  range  h  >  2  ;  in  terms  of  the  data  density 

h  correlation  length  ,  c 

p— —  t  VJ.DJ 

distance  between  data  points 

this  means  that  M  decreases  slowly  for  p  >  2  .  It  is  evident 
that  M  assumes  a  maximum,  which  equals  the  full  variance,  if 
(no  data)  and  a  minimum,  which  is  equal  to  zero,  if 


0 


30 


r  -*•  ®  (total  data  coverage)  . 

h)  In  the  range  0.5  (o> 0.5)  the  maximum  interpolation  er¬ 

ror  decreases  with  increasing  a  . 

c)  It  is  remarkable  that,  even  with  the  a  =  0.5  covariance  func¬ 
tion,  M  is  as  low  as  i  h  of  the  square  root  of  the  variance 
for  c  =  4  . 

d)  The  choice  of  the  covariance  function  becomes  extremely  criti¬ 
cal  for  ^  2  because  of  the  large  differences  in  the  interpo¬ 

lation  error  estimates.  (This  statement  is  only  valid  if  the 
data  are  free  of  noise  as  will  be  shown  later.) 

e)  The  Gaussian  covariance  function  gives  by  far  too  optimistic 
estimates . 

A  practical  example  may  illustrate  that:  assume  a  variance  of 
500  mgal2  ,  a  correlation  length  of  40  km  ,  a  reciprocal  distance 
type  covariance  function  (n  =  V2  )  ,  and  a  gravity  data  density 

of  c  =  4  (station  distance  =  10  km)  ;  on  the  basis  of  these  data 
the  maximum  interpolation  error  is  on  the  order  of  ;0.5  mgal  ;  a 
station  distance  of  20  km  leads  to  an  estimate  of  ±4  mgal. 

It  is  worth  being  mentioned  that  for  a  data  density  up  to 
-  4  the  error  estimate  depends  primarily  on  the  very  close 
neighborhood  of  the  interpolation  point;  quite  a  good  estimate  can, 
therefore,  be  obtained  by  simply  considering  the  contribution  of 
the  2  closest  data  points.  Equation  (5.1)  yields 


M2  -  *(0) 


f  ( 0 )  —  ^  (  a  ) 


(5.6) 


The  picture  of  maximum  interpolation  error  behavior  changes 
drastically  if  data  noise  is  introduced.  Let  us  first  investigate 
how  the  sampling  function  changes  if  the  data  are  superimposed  by 


'JO 


mutually  uncorrelated  errors  with  an  error  variance  t  ;  we  as¬ 
sume  that  <<•  >  (0)  .  For  a  sufficiently  stable  covariance 

matrix  C  of  exact  data,  the  incorporation  of  noise  will  change 
the  corresponding  inverse  only  liitle;  therefore,  a  linearization 
of  the  inversion  process  can  be  made  and  the  inverse  covariance 
matrix  of  the  data,  superimposed  by  noise,  is  given  by 

(C  +  cl)"1  =  C"1 (1  -  eC-i)  . 

In  terms  of  sampling  functions  o(x)  and  base  functions  ?(x)  , 

this  means  that,  with 

vfx-i)  :  =U(i-j)  +  e  5  i  _  j  ~  1  $  ( x-  j  )  ,  (5.7) 

.  is  given  by 

v(x-i)  =  v(x-i)  -  e  f  i)  (  i-  j  )  !  •»  (x-  j  )  .  (5.7)' 

Since  ,  is  a  sampling  function,  it  follows  that 
AO)  -  l  ~ 

where  '  denotes  the  diagonal  element  of  the  inverse  covariance 

matrix.  We  know  that  ;7 ‘  is  a  positive  quantity  (this  follows 

from  the  property  of  the  covariance  matrix  considered  here) ;  is 

als  positive  (a  priori  data  errcr  variance) .  Therefore, 

„(0)  •  „(0)  =  1  .  (5.3) 

It  can  be  shown  that  *  flattens  out  with  increasing 
noisy  data,  the  prediction  error  is  given  bv 


i-  or 
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nv  (x)  =  •!>  (O)  -  ~  v  (x-’<)  > (x-k)  .  (9.9) 

Introducing  (9.7)'  for  j  we  obtain 

n-(x)  ^  MO)  -  £  v  (  x-k)  ?  (x-k)  +  c  7^  :  “  ‘  (k- j  )  y  (x-j  )  <i(x-k)  ; 

k  k  j 

observing  (5.1) '  and  the  symmetry  properties  of  v  and  t,  ,  the 
prediction  error  based  on  noisy  data  is  given  by 

eo 

m2(x)  =  m:(x)  +  e  l  v2(x-k)  (5.10) 

1 2  is  non-negati ve;  therefore, 

m; (x)  >  m2 (x)  v  x  (5.11) 

as  to  be  expected.  Fig.  5.2  shows  m(;/;)  ,  dependent  on  the  data 

density,  for  various  covariance  models.  The  variance  was  put  equal 
to  1  ,  the  a  priori  data  error  variance  equal  to  0.01  . 

The  graphs  show  very  clearly  that  the  prediction  error  esti¬ 
mation  is  much  less  sensitive  with  respect  to  the  covariance  func¬ 
tion  if  noise  is  introduced  (cf.  Fig.  5.1).  In  this  particular 
case  ( e  =  1  %  of  the  variance)  there  is  practically  no  gain  in 
accuracy  achievable  if  the  data  density  o  is  greater  than  3  . 

A  practical  example  (the  same  as  before,  here,  however,  with  non¬ 
vanishing  c  =  0.01  •  varlAg}  =  5  mgal2)  may  illustrate  that:  with 
a  gravity  data  density  o  =  2  (station  distance  =  20  km)  a  pre¬ 
diction  error  just  in  between  the  data  is  estimated  as  i4.4  mgal; 
'note  that  this  value  is  only  10  %  higher  compared  to  the  =  0 
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FIG.  5.2  Estimation  of  the  prediction  error  m(V2) 

case  -  the  interpolation  error  contributes  90  %  of  the  figure)  . 

If  the  station  distance  is  reduced  to  10  km  (p  =  4)  ,  the  pre¬ 
diction  error  drops  to  ±2  mgal  and  does  not  decrease  significant¬ 
ly  with  increasing  data  density. 

Quite  a  remarkable  phenomenon  can  be  observed  in  connection 
with  least-squares  prediction  based  on  regularly  distributed  noisy 
homogeneous  data  :  under  certain,  but  certainly  not  unusual  circum¬ 
stances,  the  prediction  error  at  non-observed  points  can  be  smaller 
than  the  prediction  error  at  the  data  points.  In  the  case  discussed 
above  ( k  =  0.01  •  var{Ag}  )  ,  this  phenomenon  can  already  be  ob¬ 

served,  if  the  station  distance  is  smaller  than  12  km  (o  >  3.3)  . 
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The  reason  is  the  following:  at  the  data  points  there  is  no  pre¬ 
diction,  but  only  filtering  involved;  therefore,  the  "prediction" 
error  stems  only  from  data  noise.  In  between  the  data  interpolation 
and  error  propagation  contribute  to  the  prediction  error.  For  suf¬ 
ficiently  large  correlation  length  (or  equivalently,  for  sufficient¬ 
ly  high  data  density)  the  interpolation  error's  contribution  be¬ 
comes  insignificant  (actually,  it  goes  to  zero  if  o  -  «.  )  ;  the 
prediction  error  stems  almost  exclusively  from  the  error  propaga¬ 
tion.  If  the  errors  are  not  correlated,  the  propagated  error  is, 
indeed,  smaller  than  the  individual  error  (e.g.  error  of  the  mean 
value) .  This  is  why  such  a  curious  phenomenon  can  be  observed  in 
least-squares  prediction. 
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